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Abstract 

In  this  thesis,  the  application  of  numerical  techniques  to  electromagnetic  problems  in 
microelectronic  and  radar  imaging  systems  are  investigated.  In  particular  the  follow¬ 
ing  problems  are  studied:  (1)  Dielectric  rib  waveguide  discontinuities  are  analyzed 
with  the  Finite  Difference  Time  Domain  (FDTD)  method.  The  application  of  Beren- 
ger’s  Perfectly  Matched  Layer  to  multi-layered  dielectrics  is  analyzed  and  the  specific 
conditions  needed  to  successfully  match  the  multiple  dielectric  layers  are  determined 
and  justified.  An  FDTD  method  to  find  the  fundamental  mode’s  spatial  distribution 
is  used  to  excite  the  discontinuity  problem.  It  is  shown  that  the  computational  do¬ 
main  can  be  reduced  by  twenty  percent  over  Gaussian  excitations.  The  effects  of  rib 
waveguide  bend  discontinuities  and  the  effects  of  the  rib  geometry  to  the  bend  loss  are 
presented.  (2)  An  Impedance  Boundary  Condition  (IBC)  for  two  dimensional  FDTD 
simulations  containing  thin,  good  conductor  sheets  is  developed.  The  IBC  uses  a 
recursive  convolution  scheme  based  on  approximating  the  conductor’s  impedance  as 
a  sum  of  exponentials.  The  effects  of  FDTD  parameters  such  as  grid  size  and  time 
step  on  simulation  accuracy  are  presented.  The  IBC  is  shown  to  accurately  model 
the  conductor  loss  over  a  wide  frequency  range.  The  verification  is  performed  by 
comparing  the  quality  factors  of  rectangular  resonant  structures  determined  by  the 
FDTD  simulation  and  analytical  methods.  (3)  Phase  unwrapping  techniques  for  the 
inversion  of  terrain  height  using  Synthetic  Aperture  Radar  Interferometry  (InSAR) 
data  are  analyzed.  The  weighted  least  squares  and  branch  cut  phase  unwrapping 
techniques  are  specifically  studied.  An  optimal  branch  cut  method  and  a  hybrid  least 
squares/branch  cut  method  are  presented  and  used  to  unwrap  the  phase  of  both  sim¬ 
ulated  and  real  SAR  interferograms.  When  used  to  invert  terrain  height,  these  new 
SAR  phase  unwrapping  methods  offer  over  fifty  percent  reduction  in  root  mean  square 
(rms)  height  error  compared  to  the  straight  least  squares  method  and  over  thirty-five 
percent  reduction  in  rms  height  error  compared  to  the  weighted  least  squares  method 
based  on  coherence  data  weighting  schemes. 
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Chapter  1 
Introduction 


Numerical  modeling  of  electromagnetic  phenomena  is  a  popular  pursuit  for  many  sci¬ 
entists  and  engineers.  One  of  the  most  popular  time  domain  numerical  techniques  is 
the  Finite  Difference  Time  Domain  Method  or  FDTD.  The  FDTD  method  attempts  to 
model  the  interaction  between  electric  and  magnetic  fields  by  describing  the  universe 
at  discrete  locations  at  discrete  times.  The  physics  are  defined  by  a  set  of  difference 
equations  that  describe  the  relationships  between  field  quantities  at  the  various  loca¬ 
tions  at  different  times.  If  space  and  time  are  divided  into  small  enough  intervals  the 
difference  equations  adequately  approximate  the  electromagnetic  interactions. 

The  key  to  electromagnetic  interactions  is  the  media  and  the  FDTD  method  can 
model  most  materials  quite  accurately.  However,  computer  resources  always  constrain 
the  types  of  problems  that  can  be  simulated  within  a  reasonable  amount  of  time.  In 
fact,  FDTD’s  most  serious  limitation  is  that  it  requires  large  amounts  of  computer 
memory. 

The  computer  resources  needed  to  simulate  an  electromagnetic  problem  are  usu¬ 
ally  directly  related  to  the  smallest  wavelength  of  the  simulation.  However,  a  very 
small  wavelength  may  imply  small  structures,  and  a  small  simulation  volume  may  off¬ 
set  the  high  density  of  sample  points.  The  resource  problem  is  most  pronounced  when 
a  small  portion  of  the  simulation  contains  a  dense  material  that  greatly  shortens  the 
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wavelength.  In  this  case,  the  dense  material  drives  the  spatial  interval  and  in  turn  the 
memory  required.  Another  factor  that  greatly  influences  the  memory  requirements  is 
the  treatment  of  the  simulation  volume’s  physical  boundaries.  Modeling  open  space 
implies  no  physical  boundaries;  however,  a  simulation  volume  of  infinite  extend  is,  of 
course,  impossible.  Many  techniques  have  been  developed  to  simulate  open  space  with 
a  closed  simulation  volume;  however,  they  are  problem  specific.  It  is  the  engineer’s 
job  to  make  the  appropriate  approximations  to  simulate  a  given  problem  within  the 
computer  resource  constraints. 

Discretization  of  space,  time  and  Maxwell’s  equations  is  clearly  a  numerical  tech¬ 
nique  that  can  be  applied  to  many  types  of  electromagnetic  problems.  Other  numeri¬ 
cal  techniques  involve  the  analysis  of  the  measurements  of  electromagnetic  phenomena 
such  as  analyzing  radar  data.  Extracting  useful  information  from  radar  measurements 
in  the  presence  of  noise  is  a  continuing  topic  of  great  interest.  The  use  of  multiple 
passes  of  a  satellite  based  synthetic  aperture  radar  (SAR)  to  perform  interferometry 
is  a  very  useful  tool  to  image  the  earth’s  terrain  height.  Noise  in  the  system  and 
decorrelation  of  the  two  passes  add  greatly  to  the  height  errors  in  an  inverted  height 
image.  Reducing  these  errors  is  of  great  interest  in  the  remote  sensing  community. 

In  this  thesis,  the  application  of  numerical  techniques  to  electromagnetic  prob¬ 
lems  in  microelectronic  and  radar  imaging  systems  are  investigated.  In  particular  the 
following  problems  are  studied:  (1)  Dielectric  rib  waveguide  discontinuities  are  ana¬ 
lyzed  with  the  Finite  Difference  Time  Domain  (FDTD)  method.  Berenger’s  Perfectly 
Matched  Layer  application  to  multi-layered  dielectrics  is  analyzed  and  the  specific 
conditions  needed  to  successfully  match  the  multiple  dielectric  layers  are  determined 
and  justified.  An  FDTD  method  to  find  the  fundamental  mode’s  spatial  distribution 
is  used  to  excite  the  discontinuity  problem.  The  effects  of  this  type  of  excitation 
on  computational  domain  reduction  is  presented.  These  numerical  techniques  are 
used  to  study  the  effect  of  bend  discontinuities  in  the  rib  waveguide.  (2)  A  new  Im¬ 
pedance  Boundary  Condition  (IBC)  is  analyzed  and  developed  for  two  dimensional 
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FDTD  simulations  containing  thin,  good  conductor  sheets.  The  need  for  a  very  large 
computational  domain  is  overcome  by  the  replacement  of  the  good  conductor  with 
an  impedance  boundary  condition  that  accurately  models  the  conductor  loss  over  a 
wide  frequency  range.  The  IBC  uses  a  recursive  convolution  scheme.  The  two  di¬ 
mensional  IBC  is  verified  by  comparing  the  quality  factors  of  rectangular  resonant 
structures  determined  by  the  FDTD  simulation  and  those  calculated  with  analytical 
methods.  (3)  Phase  unwrapping  techniques  for  the  inversion  of  terrain  height  using 
Synthetic  Aperture  Radar  Interferometry  (InSAR)  data  are  analyzed.  The  weighted 
least  squares  and  branch  cut  phase  unwrapping  techniques  are  specifically  studied. 
An  optimal  branch  cut  method  and  a  hybrid  least  squares/branch  cut  method  are 
presented  and  used  to  unwrap  the  phase  of  both  simulated  and  real  SAR  interfero- 
grams.  These  new  SAR  phase  unwrapping  methods  are  compared  to  straight  least 
squares  unwrapping  and  other  weighted  least  squares  schemes. 

1.1  Technical  Discussion 

1.1.1  Analysis  of  Dielectric  Waveguide  Discontinuities 

Dielectric  waveguides  are  commonly  used  as  interconnects  in  millimeter-wave  and 
sub-millimeter  wave  integrated  circuit  technologies  [75]- [100].  As  frequencies  increase, 
conductor  losses  diminish  the  utility  of  microstrip  and  coplanar  interconnect  struc¬ 
tures.  Furthermore,  to  maintain  single  mode  operation  the  structure  size  must  de¬ 
crease  which  further  increases  the  loss  through  these  metal  guides.  Dielectric  waveg¬ 
uides  avoid  these  losses;  however,  when  these  interconnects  contain  discontinuities 
such  as  bends,  these  discontinuities  introduce  loss  through  the  waveguide.  Since  the 
propagating  modes  of  this  type  of  structure  are  characterized  by  very  complex  field 
distributions  that  do  not  have  exact  analytic  solutions,  a  numerical  approach  is  ap¬ 
propriate  to  investigate  these  structures.  In  this  thesis  an  FDTD  implementation  is 
used  to  analyze  the  effects  of  waveguide  discontinuities  on  the  fundamental  mode’s 
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propagation  through  the  waveguide. 

The  main  difficulty  of  using  the  FDTD  method  on  this  type  of  problem  is  computer 
resources.  With  multiple  dielectric  media,  the  most  dense  material  determines  the 
spatial  interval  of  the  discretization.  And,  the  open  nature  commonly  found  in  many 
dielectric  guiding  structures  creates  the  need  for  special  treatment  on  the  edges  of 
the  simulation  volume. 

The  first  analysis  of  dielectric  waveguides  focused  on  finding  dispersion  curves 
and  field  distributions  in  two  dimensional  cross  sections  of  the  guides  using  frequency 
domain  techniques  [75,  76,  77,  78].  More  recent  work  has  concentrated  on  using 
integral  equations  to  solve  the  two-dimensional  problem  [82,  99],  while  others  have 
used  frequency  methods  like  Beam  Propagation  Method  or  the  discretization  of  the 
scalar  wave  equation  [81,  83].  Much  of  the  dielectric  waveguide  discontinuities  work 
has  been  limited  to  two  dimensional  dielectric  slab  waveguides  as  in  [94]. 

Three  dimensional  single  frequency  analysis  has  been  done  on  structures  such 
as  directional  couplers  [80],  tapered  rib  waveguides  [89],  Y-junctions  [92],  and  step 
discontinuities  [84]. 

Multi-frequency  analysis  of  three  dimensional  structures  using  FDTD  has  been 
done  on  the  transitions  from  rectangular  waveguides  and  microstrip  to  dielectric 
waveguides  [85,  96,  100],  but  little  has  been  done  with  dielectric  bends.  One  work 
reports  a  full  three  dimensional  look  at  a  bend  in  a  rectangular  dielectric  guide  limited 
to  a  single  frequency  [93]. 

Figure  1-1  is  a  graphical  representation  of  a  dielectric  rib  waveguide  and  the 
essential  elements  of  a  FDTD  simulation.  Unless  bounded  with  a  perfect  conductor, 
the  walls  of  the  simulation  volume  must  contain  some  absorbing  boundary  condition. 
The  problem  is  excited  by  an  excitation  plane  and  the  power  is  measured  at  two 
reference  planes. 

When  using  FDTD  to  study  any  waveguide,  depending  on  the  excitation  used,  it 
can  take  some  distance  for  the  fundamental  mode  to  develop  as  the  evanescent  modes 
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Figure  1-1:  Dielectric  rib  waveguide  configuration  in  a  3D  FDTD  simulation. 
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decay  away  and  radiating  modes  leave  the  simulation  volume.  In  order  to  more 
efficiently  launch  the  fundamental  mode,  in  this  work  the  mode’s  spatial  distribution 
is  determined  up  front.  A  two  dimensional  FDTD  method  is  used  to  construct  a  mode 
template  [131,  133].  The  three  dimensional  FDTD  code  is  collapsed  in  the  direction 
of  propagation  by  replacing  the  spatial  difference  equation  with  an  equation  based 
on  an  assumed  propagation  constant  and  the  corresponding  phase  difference  between 
adjacent  spatial  grid  points.  A  two  dimensional  simulation  of  the  waveguide  cross 
section  produces  the  corresponding  temporal  frequencies  of  the  propagating  modes. 
The  lowest  frequency  is  the  fundamental  mode  and  the  spatial  distribution  is  found 
by  performing  a  Fourier  transform  at  each  space  point  on  the  two  dimensional  grid 
at  the  fundamental  mode  frequency.  With  this  source  condition,  the  computational 
domain  can  be  reduced  from  the  commonly  used  spatial  Gaussian  source  condition 
by  allowing  shorter  distances  between  excitation  and  the  discontinuity  and  thus  the 
waveguide  can  be  studied  with  fewer  computer  resources. 

Unlike  metallic  waveguides,  the  dielectric  rib  waveguide  is  an  open  structure  and 
analyzing  them  with  the  Finite  Difference  Time  Domain  method  requires  absorbing 
boundary  conditions  (ABC)  that  simulate  open  space.  In  this  work  the  new  Perfectly 
Matched  Layer  [46]  is  used.  A  careful  examination  of  its  application  to  multi-layer 
dielectrics  is  given  and  the  results  are  used  to  implement  the  ABC  for  the  study  of 
waveguide  discontinuities. 

1.1.2  Modeling  of  Thin  Finite  Conductivity  Sheets 

In  FDTD  simulations  involving  highly  conductive  materials,  such  as  the  metal  case  of 
a  computer  system,  the  tangential  electric  fields  are  typically  set  to  zero  on  the  surface 
of  the  material.  This  Perfect  Electric  Conductor  (PEC)  assumption  ignores  any  loss 
associated  with  the  less  than  infinite  conductivity.  The  errors  that  result  from  this 
approximation  are  considered  against  the  large  cost  of  discretizing  the  lossy  material. 
Since  the  wavelength  of  a  highly  conductive  material  is  very  small,  the  simulation 
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Figure  1-2:  Replacement  of  thin  conductor  with  Impedance  Boundary  Condition 


volume  must  be  divided  into  a  very  fine  grid  in  order  to  capture  the  loss  mechanism 
within  the  conductor.  When  a  uniform  grid  is  used  (called  the  brute  force  method) 
computer  resources  are  wasted  by  having  to  finely  divide  the  other  less  dense  mate¬ 
rials  {e.g.  free  space).  In  these  materials  the  fine  grid  is  not  necessary  to  capture  the 
physics  of  the  problem.  The  computational  size  can  be  reduced  with  a  sub-gridding 
scheme  where  the  conducting  material  uses  a  much  smaller  grid  than  the  rest  of  the 
computational  domain  [9,  10,  12].  With  this  technique  care  must  be  given  to  reduce 
reflections  caused  by  the  change  in  lattice  grid  size.  For  very  good  conductors;  how¬ 
ever,  the  grid  size  must  be  so  small  that  even  sub-gridding  is  not  a  viable  option. 
To  overcome  the  resource  problem,  the  surface  of  the  highly  conductive  material  can 
be  replaced  with  an  Impedance  Boundary  Condition  (IBC)  [126]-[130].  Figure  1-2 
represents  how  the  conductive  material  is  replaced  with  a  boundary  condition  that 
incorporates  the  physics  of  the  conductor  layer.  An  IBC  is  only  appropriate  when  the 
simulation  volume  of  interest  is  on  one  side  of  the  the  conductive  material.  However, 
IBCs  have  the  added  complication  that  they  are  usually  frequency  dependent  and  are 
not  directly  applicable  to  the  the  standard  frequency  independent  FDTD  equations. 
In  this  case,  the  FDTD  equations  must  be  modified  to  incorporate  the  dispersive 
nature  of  the  surface  [109]  -[120].  Typical  frequency  domain  equations  relating  the 
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electric  fields  and  magnetic  fields  become  convolution  equations  in  time.  Since  a 
convolution  has  a  large  memory  overhead,  if  the  convolution  integral  can  be  approxi¬ 
mated  by  a  sum  of  exponentials,  then  recursive  convolution  can  be  implemented  and 
the  memory  requirements  reduced. 

Similar  to  an  IBC,  another  approach  [117]  uses  a  synthetic  conductivity  and  the 
normal  FDTD  difference  equations  for  lossy  media.  A  synthetic  conductivity  is  de¬ 
rived  by  comparing  the  numerical  impedance  of  the  difference  equations  and  the 
actual  impedance  of  a  good  conductor  at  specified  frequency.  In  this  way  the  derived 
synthetic  conductivity  is  inserted  into  the  FDTD  equations  at  the  boundary;  thus 
being  a  surface  boundary  condition  similar  to  an  IBC.  The  advantage  to  this  method 
is  that  no  new  equations  are  needed;  however,  it’s  major  disadvantage  is  that  it  is 
only  appropriate  at  a  single  frequency. 

An  Impedance  Boundary  Condition  that  incorporates  frequency  dependence  was 
developed  for  thin  dielectric  coatings  over  PEC  surfaces  [127].  The  IBC  formulation 
starts  with  the  analytically  derived  expression  for  the  impedance  in  the  frequency 
domain.  The  expression  is  expanded  in  a  Taylor  series  about  the  wave  number.  Next 
the  frequency  domain  is  transformed  to  the  time  domain  by  replacing  -iu  with  d/dt. 
The  resulting  partial  differential  equation  is  fourth  order  in  time  and  first  order  in 
space.  The  difference  equations  do  model  the  frequency  characteristics  fairly  well; 
however,  the  equations  are  quite  complicated  and  rely  on  many  past  values  of  the 
surface  electric  field.  Furthermore,  the  condition  is  limited  to  normal  incidence;  thus, 
limited  to  ID  structures.  This  work  is  extended  to  include  non-normal  incidence 
[129];  however,  the  modified  IBC  is  only  good  for  a  single  incidence  angle. 

Like  the  IBC  of  [127],  Tassoudji  [128]  starts  with  a  frequency  domain  expression 
for  the  impedance  of  one  dimensional  structures.  Here  the  structure  is  a  thin  sheet 
of  highly  conducting  material.  The  frequency  domain  equation  is  transformed  to  the 
Laplace  domain  and  the  subsequent  expression  is  approximated  by  a  sum  of  first- 
order  rational  functions.  This  approximation  is  critical  to  the  IBC’s  practicality  to 
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Figure  1-3:  SAR  interferogram 

FDTD.  The  first-order  rational  functions  are  exponentials  in  time  and  thus  a  recursive 
convolution  scheme  is  implemented.  The  benefit  of  this  scheme  is  the  wide  range  of 
frequencies  for  which  it  is  valid. 

In  this  thesis,  the  IBC  for  thin  finite  conducting  sheets  of  [128]  is  examined  and 
extended  to  two  dimensions.  Guidelines  for  the  number  of  expansion  terms  needed  are 
derived  in  order  to  give  users  an  a  priori  estimate  of  the  resources  required  to  perform 
a  given  simulation.  The  new  2D  IBC  is  validated  through  the  measurements  of  the 
quality  factor  of  rectangular  resonant  cavities  at  different  frequencies  and  conductivity 
values. 


1.1.3  Phase  Unwrapping 

Airborne  and  spaceborne  Synthetic  Aperture  Radar  (SAR)  platforms  have  been  used 
for  many  years  to  study  the  earth’s  surface  [144]-[159].  When  two  radars  on  a  single 
platform  or  two  passes  of  a  single  radar  map  the  same  area,  an  interferogram  can  be 
produced  from  the  difference  in  phase  measured  by  each  radar  or  pass.  An  interfero¬ 
gram  is  a  pictorial  representation  of  the  phase  differences  measured  at  each  pixel  as 
shown  in  Figure  1-3 

Since  the  measured  phase  differences  lie  between  — tt  and  tt,  the  phase  is  said  to 
be  wrapped.  A  SAR  interferogram  contains  fringes.  These  fringes  are  the  locations 
on  the  interferogram  where  a  27r  discontinuity  exists.  The  interferogram  resembles  a 
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Figure  1-4:  Digital  Elevation  Model  (DEM)  of  Alaskan  mountain. 

topographical  contour  map  where  a  line  of  constant  elevation  corresponds  to  a  fringe. 
When  no  noise  is  present,  the  fringes  can  easily  be  located  and  the  data  adjusted  by 
adding  multiples  of  2tt  to  produce  an  unwrapped  phase  image.  However,  real  data 
contain  noise  and  phase  unwrapping  can  be  a  complicated  process. 

Successful  phase  unwrapping  is  the  key  to  the  extraction  of  DEM  (digital  elevation 
model)  from  an  interferometric  SAR  phase  image.  Figure  1-4  is  the  DEM  associated 
with  the  terrain  that  produced  the  ERS-1  SAR  interferogram  of  Figure  1-3. 

Although  some  new  phase  unwrapping  techniques  have  been  introduced  recently 
based  on  the  principle  of  maximum  entropy  [164]  and  multiresolution  [165]  that  may 
lead  to  better  phase  unwrapping  algorithms,  there  are  really  two  basic  approaches  to 
unwrapping  phase  data.  The  first  is  based  on  finding  an  unwrapped  solution  such 
that  the  solution’s  first-order  partial  derivatives  in  the  x  and  y  directions  match  (or 
closely  match)  the  wrapped  first-order  partial  derivatives  or  gradients  of  the  phase 
data.  Typically  noise  is  handled  by  unwrapping  the  best  data  first  in  local  schemes 
or  weighting  the  data  in  global  schemes.  The  second  approach  integrates  along  the 
data  and  adds  or  subtracts  2n  when  a  fringe  is  crossed.  The  noise  is  considered  by 
searching  for  phase  inconsistencies,  in  the  form  of  residues.  The  residues  are  connected 
to  form  branch  cuts  and  phase  is  unwrapped  by  adjusting  the  integration  path  or  by 
modifying  the  fringe  information. 
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There  are  local  and  global  techniques  to  phase  unwrapping  algorithms  based  on 
the  gradients  of  the  measured  data.  The  local  algorithms  usually  locate  one  or  more 
areas  on  the  image  that  are  considered  good  data.  Finding  these  points  may  use  noise 
floor  data  and  signal-to-noise  ratio  [168]  or  use  coherence  data  [166], [178].  Unwrap¬ 
ping  usually  begins  with  these  good  areas.  The  algorithms  move  from  pixel  to  pixel 
and  add  or  subtract  27r  based  on  some  criteria.  One  method  follows  the  least-gradient 
path  with  the  assumption  that  the  smallest  gradient  points  to  the  best  data  [180], [178] 
and  the  adjustments  are  made  to  match  the  solution  gradient  to  the  wrapped  mea¬ 
sured  gradient.  Other  methods  use  more  neighbors  to  decide  the  value  to  add  to  the 
unwrapped  pixel  [161], [169].  Here  all  the  gradients  of  the  adjacent  unwrapped  pixels 
are  examined  and  the  phase  of  the  unwrapped  pixel  is  predicted  based  on  these  gra¬ 
dients.  Then  an  integer  number  of  27r  is  added  to  bring  the  unwrapped  pixel  closest 
to  the  predicted  phase.  In  this  way  it  is  possible  to  have  phase  differences  larger  than 
TT,  so  this  method  can  accommodate  real  discontinuities  due  to  terrain  features.  All 
of  these  local  techniques  have  the  added  complication  of  growing  separate  unwrapped 
regions  that  must  be  joined  to  produce  the  final  product.  Since  the  best  pixels  are 
unwrapped  first,  the  likelihood  of  errors  propagating  through  the  image  is  reduced. 

Global  gradient  algorithms  are  based  on  least  squares  and  weighted  least  squares 
methods  [170].  In  this  way,  the  technique  attempts  to  find  a  solution  that  minimizes 
the  differences  between  all  of  the  solution’s  gradients  and  the  wrapped  data’s  gradi¬ 
ents.  The  least  squares  approach  is  very  desirable  because  of  the  speed  at  which  all 
the  data  can  be  unwrapped;  but,  the  method  does  not  treat  noisy  data  very  well  [170]. 
The  weighted  least  squares  method  allows  the  user  to  favor  the  good  data  by  applying 
a  set  of  weights  to  the  data  based  on  some  knowledge  of  the  data’s  noise  content. 
These  methods  were  applied  to  a  single  SAR  interferogram  [160];  however,  no  details 
of  the  weighting  criteria  were  presented.  [162]  used  a  weighting  least  squares  on  a  sim¬ 
ulated  interferogram  with  a  shear,  where  the  shear  discontinuity  was  masked  out  and 
a  new  multi-grid  iteration  scheme  was  used  to  solve  the  weighted  partial  differential 
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equation.  [171]  applied  a  weighted  least  squares  technique  to  a  speckle  interferogram 
using  a  weighting  scheme  based  on  masking  out  the  phase  inconsistencies  or  residues. 

The  branch  cut  method  also  has  local  and  global  approaches.  So  far,  only  local 
approaches  have  been  reported  in  SAR  interferometry.  The  branch  cut  method  and 
SAR  interferometry  phase  unwrapping  was  first  reported  in  [173].  Applied  to  the 
interferogram  derived  from  two  passes  of  Seasat,  the  residue  connections  were  based 
on  a  nearest  neighbor  approach.  This  approach  works  well  with  a  low  density  of 
residues;  however,  it  breaks  down  with  a  high  density  of  residues  [173].  Once  two 
residues  are  connected  they  are  removed  and  no  longer  considered  for  connection  to 
any  other  residues.  This  method  is  very  likely  to  leave  an  uncompensated  residue. 
The  phase  is  then  unwrapped  by  integrating  along  a  path  that  never  crosses  a  branch 
cut.  This  approach  was  briefly  reported  in  [160]  as  part  of  an  overview  of  SAR  phase 
unwrapping  techniques.  The  reported  disadvantage  to  the  branch  cut  method  is  the 
propagation  of  global  errors  from  uncompensated  residues.  A  modification  to  the 
basic  branch  cut  method  in  SAR  interferometry  involves  connecting  and  removing 
residue  pairs  that  are  only  1  or  2  pixels  apart  [163].  The  remaining  residues  are 
handled  separately  and  considered  part  of  real  discontinuities  that  exist  as  a  result  of 
terrain  features. 

The  branch  cut  method  has  also  been  applied  to  speckle  interferometry,  used 
to  measure  very  small  surface  deformations  on  structures.  In  this  application,  a 
nearest  neighbor  (local)  connection  algorithm  was  used  and  reported  in  [174].  The 
first  global  branch  cut  method  is  used  to  unwrap  a  speckle  interferogram.  All  residues 
are  considered  before  making  any  branch  cuts.  In  this  way,  the  algorithm  is  based  on 
minimizing  the  total  branch  cut  length[175]. 

Finding  efficient,  accurate  and  automatic  phase  unwrapping  algorithms  is  an  inter¬ 
esting  topic  and  is  becoming  more  important  as  more  processing  occurs  on  the  SAR 
platform.  In  this  thesis,  the  least  squares  global  approaches  are  investigated  and 
applied  to  SAR  interferometry.  Specifically,  a  hybrid  method  that  uses  a  weighting 
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mask  based  on  branch  cut  residues  is  examined  and  compared  to  weighting  schemes 
based  on  coherence  data.  Also,  optimum  branch  cut  algorithms  are  applied  to  SAR 
interferograms  and  the  techniques  are  compared  to  the  least  squares  method.  Data 
from  both  simulated  and  real  SAR  interferograms  are  used. 


1.2  Description  of  Thesis 

This  thesis  is  divided  into  six  chapters.  Chapters  1  and  6  are  the  introduction  and 
conclusion  and  require  no  description. 

Chapter  2  provides  a  detailed  look  at  the  Finite  Difference  Time  Domain  Method 
as  a  numerical  technique  to  solving  electromagnetic  problems.  The  discretization 
of  Maxwell’s  equations  and  the  treatment  of  both  lossless  and  lossy  materials  are 
described.  Special  attention  is  given  to  Absorbing  Boundary  Conditions  and  specifi¬ 
cally  the  relatively  new  perfectly  matched  layer  and  the  it’s  application  to  multilayer 
dielectric  media  in  preparation  for  it’s  use  in  Chapter  3. 

Chapter  3  is  dedicated  to  the  use  of  the  numerical  FDTD  technique  to  study 
dielectric  rib  waveguides.  The  perfectly  matched  layer  of  Chapter  2  is  used  as  the 
absorbing  boundary.  The  application  of  a  mode  template  to  efficiently  launch  the 
guide’s  fundamental  mode  is  described.  This  numerical  technique  is  compared  to 
analytical  methods  for  deriving  dispersion  characteristics  of  dielectric  waveguides. 
A  numerical  study  of  the  benefit  of  this  technique  is  conducted.  Afterwards,  these 
methods  are  used  to  analyze  various  bend  discontinuities  found  in  dielectric  waveguide 
structures. 

In  Chapter  4  a  new  Impedance  Boundary  Condition  is  extended  for  use  in  two 
dimensional  FDTD  simulations.  The  IBC  is  described  in  great  detail  and  studied  in 
one  dimensional  simulations.  The  influence  of  FDTD  parameters  and  the  number  of 
IBC  expansion  terms  on  the  accuracy  of  the  IBC  based  on  the  reflection  of  a  wave  from 
a  thin  finite  conducting  sheet  is  presented.  The  IBC  is  extended  to  two  dimensions 
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and  verified  through  the  study  of  resonant  structures.  The  quality  factor  of  many 
rectangular  resonant  structures  is  found  using  FDTD  and  compared  to  analytical 
results. 

Numerical  techniques  of  phase  unwrapping  of  synthetic  aperture  radar  interfero- 
grams  are  investigated  in  Chapter  5.  The  foundations  of  global  phase  unwrapping  are 
described,  specifically  least  squares  and  optimal  branch  cut  unwrapping.  New  mod¬ 
ifications  to  these  techniques  are  applied  to  both  simulated  and  real  interferograms. 
After  the  application  of  these  techniques  to  the  real  SAR  data,  the  phase  is  used  to 
invert  terrain  height  and  the  height  is  compared  to  ground  truth. 

There  are  three  appendices  included  in  this  thesis.  Appendix  A  lists  the  first 
and  second  order  Mur  absorbing  boundary  condition  FDTD  equations.  Appendix  B 
describes  the  methods  used  to  determine  the  coefficients  used  in  the  PML  FDTD 
equations.  Appendix  C  gives  a  detailed  review  of  the  transportation  problem  found 
in  linear  programming  that  was  used  to  build  an  optimum  branch  cut  connection 
algorithm  in  Chapter  5. 


Chapter  2 


Finite  Difference  Time  Domain 
Method 


2.1  Introduction 


The  Finite  Difference  Time  Domain  Method  is  one  of  the  most  popular  numerical 
methods  for  solving  electromagnetic  problems  in  the  time  domain.  K.  S.  Yee  [1]  pro¬ 
posed  a  discretizing  scheme  for  Maxwell’s  equations  that  staggers  the  electric  and 
magnetic  fields  in  space  and  time  in  such  a  way  that  the  difference  equations  are 
second  order  accurate.  But,  more  importantly,  it  is  possible  to  analyze  very  com¬ 
plicated  structures  over  multiple  frequencies  with  a  single  simulation.  FDTD’s  basic 
limitation  is  its  need  for  large  computer  resources  to  represent  the  simulation  volume. 
The  impact  of  this  limitation  has  been  lessened  with  the  recent  advances  in  computer 
technology.  The  FDTD  method  has  been  widely  used  to  solve  electromagnetic  prob¬ 
lems  [2]-[28].  The  purpose  of  this  chapter  is  to  describe  the  FDTD  method  used  in 
this  thesis. 


35 


36 


CHAPTER  2.  FINITE  DIFFERENCE  TIME  DOMAIN  METHOD 


2.2  Maxwell’s  Equations  in  Cartestian  Coordinates 


Maxwell’s  curl  equations  form  the  basis  of  the  FDTD  difference  equations.  Starting 
with 

V  X  E{f,  t)  =  t),  (2.1) 

and 

d 

V  X  H{r,t)  = -^D{r,t)  + (2.2) 

with  the  constitutive  relations 

D  =  eE,  (2.3) 

and 

B  =  f^H,  (2.4) 

the  first  step  is  to  separate  the  equations  into  their  individual  components.  The 
electric  field  component  equations  are 


and 


^  ,  IP  ^  ^  ZJ  ^  rr  J 

-^H  -  T 

dt^  ~ 


^  eE..  =  S-H..  -  ^H^  -  J„ 


dt 


dx  dy 


while  the  magnetic  field  equations  are 


(2.5) 

(2.6) 

(2.7) 


l„r,  -Af 

"  dz^’  dy^‘’ 


A„H  -Af  -Af 

dA^’'~ax^‘  dz^" 
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(2.9) 
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and 


“  dy^’  dx^’’ 


(2.10) 


2.3  Difference  Equations  in  Free  Space 


2.3.1  Central  Difference  Approximation 

The  FDTD  difference  equations  are  constructed  from  a  central  difference  approxima¬ 
tion  of  the  partial  derivatives  in  equations  (2.5)  through  (2.10). 

Given  a  function  f(x),  the  continuous  variable  x  is  divided  into  discrete  points 
separated  by  A  such  that  x  =  iA  where  i  is  any  integer.  The  partial  derivative  may 
be  approximated  by 


dx 


f{^) 


x=iA 


f{(i  +  l)A)-f{{i-^,)A) 
A 


(2.11) 


This  central  difference  approximation  is  second-order  accurate  (i.e.  the  error  ~  A^). 
Thus,  theoretically,  the  error  can  be  made  as  small  as  desired  for  well  behaved  func¬ 
tions  if  the  grid  size.  A,  is  made  small  enough. 


2.3.2  Yee  Grid  and  FDTD  Nomenclature 

Yee  [1]  introduced  a  spatial  arrangement  of  the  six  field  quantities  shown  in  Figure 
2-1,  that  places  the  electric  field  unknowns  along  the  edges  of  a  discretized  grid  and 
the  magnetic  field  unknowns  on  the  faces  of  the  grid.  The  electric  and  magnetic  fields 
are  staggered  both  in  time  and  space  by  ^  and  and  respectively.  This 

staggered  approach  allows  the  difference  equations  to  be  local  in  both  time  and  space. 
By  local  we  mean  that  an  unknown  is  a  function  of  itself  and  its  nearest  neighbors 
(local  in  space)  at  the  most  resent  time  step  (local  in  time)  .  Thus  only  one  set  of 
unknowns  at  each  location  must  be  stored  in  memory  at  a  given  time.  The  set  of  all 
discrete  cells  is  called  the  computational  domain. 


2.3.  DIFFERENCE  EQUATIONS  IN  FREE  SPACE 


39 


The  continuous  coordinates  (x,  y,  z,  t)  become  discrete  coordinates  (i Ax,  jAy, 
kAz,  nAt)  or  simply  (i,  j,  k,  n)  where  i,  j,  and  k  represent  the  space  coordinate  indices 
and  n  represents  the  time  coordinate  index  and  Ax,  Ay,  Az  and  At  are  the  spatial 
and  temporal  intervals  between  sampled  points.  These  intervals  are  referred  to  as  the 
grid  size  (spatial)  and  time  step  (temporal).  When  uniform  gridding  is  used,  the  grid 
spacings  are  all  set  to  a  single  value,  A.  The  most  common  FDTD  nomenclature  uses 
a  subscript  to  represent  the  field  component  direction,  a  superscript  to  represent  the 
time  index,  and  three  arguments  to  represent  the  space  indices. 

For  example,  E^^^{i,j,k  +  5)  represents  the  z  component  of  the  electric  field 
sampled  at  t  =  (n  +  1)  At  seconds  and  sampled  in  space  at  x  =  iAx,  y  —  jAy  ,  and 
2:  =  (^  +  1)  Az  meters. 

2.3.3  Difference  Equations 

Figure  2-2  shows  the  spatial  arrangement  of  the  individual  field  components  used  to 
construct  the  difference  equations.  A  field  equation  includes  its  past  value  plus  a 
linear  combination  of  the  four  nearest  fields  that  lie  in  the  plane  perpendicular  to  the 
the  field  component.  The  free  space  difference  equations  are  depicted  in  equations 
2.12-2.17. 


At 


i  j  +  h  k) 


•^+0  - 


toAy 


At 

eoAz 


(i  + 


5  *:  +  5)  -  + 


(2.12) 


(2,13) 
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Figure  2-2:  Graphical  representations  used  to  construct  FDTD  difference  equations. 
Solid  lines  represent  integer  grid  lines. 
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At 


eoAz 


At 

eoAx 


^Hx  {i,j  +  ^,k  +  ^)  —  Hx  ‘^{i,j  +  ^,k  2)^ 


v- 
At 


-7^  k  +  ^)-  Hy^^{i-  \,j,  k  + 1)^ 

'1^  (^Hx^^iiJ  +  i^  +  ^)  -  HT^{i,j  -  i  A;  +  5)^ 


(2.14) 


H^hiJ  +  I,  +  i)  =  nrkiJ  +  i  ^  +  5)  (2-15) 

+7^z  +  hk  +  l)-  E^{i,j  +  i,  k)) 

-iSk  j  +  I,  k  +  ^)-  E^{i,  j,  k  +  1)) 


Hy^^{i  +  hj,k  +  ^)  =  Hy  ^(t  +  |,7,fc  +  |)  (2.16) 

(^x(*  +  i i,  k+1)-  E-{i  +  1,  j,  k  -  1)) 

(^r(* + i,i,  ^ + 5)  -  ^ + i)) 


Hz  {i  +  ^,j  +  ^,k)  =  Hz  ^(*  +  5)i  +  |)^) 
+  2>  ■?  +  +  5’ ■?>  ^)) 

+  i  +  i  ^)) 


(2.17) 
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As  shown,  the  temporal  discretization  places  the  electric  and  magnetic  fields  half  a 
time  step  away  from  each  other.  With  this  configuration,  the  method  first  calculates 
the  electric  field  (2.12-2.14),  then  advances  the  time  by  half  a  time  step  and  calculates 
the  magnetic  field  (2.15-2.17).  This  technique  is  called  leap  frogging. 

2.4  Difference  Equations  in  Lossy  Dielectric  Me¬ 
dia 

In  this  section,  we  describe  the  inclusion  of  lossy  media  and  the  scattered  versus  total 
field  FDTD  formulation. 

2.4.1  M2ixwell’s  Equations  in  Lossy  Dielectric  Media 

The  framework  for  FDTD  equations  and  algorithm  used  in  this  thesis  were  based 
on  those  in  Kunz  and  Luebber  [14]  and  modified  where  necessary.  We  begin  with 
Maxwell’s  curl  equations  for  source  free  lossy  medium: 

d  -  1  -  - 

—H  =  — VxE-—H  (2.18) 

at  p  p 

and 

1  6 

XH  +  —E.  (2.19) 

at  €  e 

When  equations  (2.18)  and  (2.19)  are  discretized,  we  have  a  total  field  formulation; 
however,  it  is  more  convenient  to  use  a  scattered  field  formulation.  The  main  benefit 
of  a  scattered  field  formulation  is  that  an  incident  field  can  be  analytically  specified 
throughout  the  computational  domain.  In  this  way  a  plane  wave  can  be  specified 
without  needing  to  provide  a  source  sufficiently  far  away  to  generate  a  plane  wave 
at  the  scatterer.  This  reduces  the  size  of  the  computational  domain  and  provides  a 
convenient  way  to  insert  plane  waves  radiating  at  different  angles.  Also,  it  is  more 
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general  than  a  total  field  formulation  i.  e.  a  scattered  field  code  becomes  a  total  field 
code  when  the  incident  field  is  set  to  zero. 

The  basic  assumption  is  that  Maxwell’s  equations  are  linear  in  the  simulation  vol¬ 
ume.  Thus,  outside  any  scatterers  both  the  incident  and  total  fields  satisfy  Maxwell’s 
equations  for  free  space  and  inside  any  scatterer  the  total  field  satisfies  Maxwell’s 
equations  for  the  lossy  medium.  Of  course,  the  linearity  assumption  does  not  hold 
for  all  media  but  does  apply  to  the  media  used  in  this  thesis.  So  we  have 


E*  =  E^  +  E\ 


(2.20) 


and 


(2.21) 


where  the  superscripts  t,  i,  and  s  stand  for  total,  incident,  and  scattered  respectively. 
The  incident  fields  satisfy  the  free  space  equations: 


V  X  E'  = 


(2.22) 


and 

VxW  =  (2.23) 

and  the  total  fields  satisfy  equations  (2.18)  and  (2.19).  Substitution  of  (2.20)  and 
(2.21)  into  (2.18)  and  (2.19)  and  subtracting  the  incident  field  equations  (2.22)  and 
(2.23)  we  find  the  scattered  field  equations  for  inside  the  scatters. 


V  X  E"  =  - 

L/  V 


[ft  -  m)^H‘  +  <rw 


(2.24) 


and 


VxH^  =  + 


{e-eo)^^E^  +  a^E^ 


(2.25) 


Outside  the  scatter,  since  the  total  fields  also  satisfy  (2.22)  and  (2.23),  substitution 


44 


CHAPTER  2.  FINITE  DIFFERENCE  TIME  DOMAIN  METHOD 


and  subtraction  of  the  incident  fields  yield 


V  xE^  =  -Po-^H^  (2.26) 

and 

V  xH^  =  eo-^E^.  (2.27) 

Equations  (2.24), (2.25), (2.26),  and  (2.27)  are  the  foundation  for  the  FDTD  difference 
equations.  Since  equations  (2.26)  and  (2.27)  are  just  special  cases  of  (2.24)  and  (2.25) 
with  fx  —  fj,Q,  e  —  eo,  cr  =0  and  cr^  =  0  only  (2.24)  and  (2.25)  are  needed  to  generate 
the  FDTD  equations.  The  equations  for  each  field  component  are 


dC 


e~El  +  a^El  = 


\d  rr.  d  1 

d 

dy  ^  dz  \ 

+ 

d  d  ,,, 

d  ■  .1 

+ 

_ 

dt 


e-Et  +  a^El  = 


d  d 


+ 


o 


JLe*  _  2.e‘ 


d 


d  .1 

1 

'  |co 

1 

- 1 

and 


(2.28) 

(2.29) 

(2.30) 

(2.31) 

(2.32) 


H^Hl  +  a^Hl  =  - 


5  d 

ipS  ^  pis 

dx  y  dy  ^ 

— 

{^x-^x,)-HiTa^Hi 


(2.33) 


For  completeness,  the  most  general  equations  are  shown  above.  Because  in  this  work 
no  magnetic  materials  are  used,  only  equations  (2.28)-(2.30)  need  to  be  discretized. 
The  magnetic  fields  can  be  updated  with  the  free  space  difference  equations  (2.15)- 
(2.17). 
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2.4.2  Difference  Equations 

The  electric  field  difference  equations  for  lossy  dielectric  media  are  shown  below. 
The  equations  are  constructed  with  computer  coding  in  mind.  The  coefficients  Ci  - 
Ci  are  calculated  up  front  and  stored  in  computer  memory  and  called  during  the 
simulation.  Also,  the  coefficients  are  a  function  of  the  medium,  m.  One  of  the  tasks 
at  the  beginning  of  the  FDTD  simulation  is  to  build  the  structure  under  test  by 
specifying  which  materials  are  at  each  cell  location;  such  that  m  is  a  function  of  the 
field  component  and  location.  The  electric  field  equations  are 

k)  =  Ct{m)E^^^{i  +  |,  j,  k)  (2.34) 

-  CI{m)Ei^{i  +  lj,k) 

-  CI{m)der{Ei^{i  +  iLk)} 

-  Cf(m)  + |)  -  5), 

E;'‘*'(iJ  +  lk)=Cl{m)E;’'{iJ  +  lk)  (2.35) 

-  C’,(m)E;’‘{i,j+lk) 

-  Ct{m)der{Ei’'{i,j  +  lk)} 

+  C'r(m)  Hr^^{i,j  +  lk  +  k)-Hr^^i,j  +  l,k-\) 

and 


Er*\i  +  \,3.  k)  =  C!(m)El’'{i  +  1 J,  k) 


(2.36) 
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-  C^,{m)Ei-{i  +  ^,j,k) 


-  CI{m)der{Ei-{i  +  lj,k)} 


+  Cr{m) 
—  C4^{m) 


nr^-^ 


{i  +  lJ,k  +  l)-Hr^'^ 


{i-k,j,k  +  l) 

{i,j  -  +  |) 


where 


Ct{m)  = 


C|(m)  = 


e 

e  +  (T®Ai  ’ 

cr^At 
£  +  cr®Ai’ 


C|(m)  = 


(e  -  fo}At 
e  +  <7®  At  ’ 


(2.37) 

(2.38) 

(2.39) 


Crim) 


At 

Aa{e  +  cr®  At) 


x,y,z, 


and  der{}  represents  a  discretized  analytic  expression  for  the  time  derivative  of  the 
analytically  specified  incident  field. 


2.5  Treatment  of  Dielectric  Interfaces 

Figure  2-3  represents  the  placement  of  the  electric  and  magnetic  fields  in  a  Yee  grid 
cell.  A  single  cell  represented  by  the  indices  {i,j,k)  has  its  lower  left  front  corner 
positioned  at  {iAx,jAy,kAz).  Although  a  cell  only  has  one  each  of  the  six  field 
quantities  associated  with  it,  a  cell  has  12  electric  and  6  magnetic  field  quantities 
adjacent  to  it  as  Figure  2-3  shows.  This  configuration  is  key  to  the  construction  of 
multiple  dielectric  media  structures  within  the  simulation  volume. 

When  more  then  one  material  is  used  in  the  simulation  volume  permittivity  be¬ 
comes  a  function  of  position.  In  general,  the  permeability  is  also  a  function  of  position, 
but  throughout  this  thesis  only  non-magnetic  materials  are  used.  Although  there  are 
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Figure  2-3:  Individual  field  locations  on  the  Yee  grid. 
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many  possible  ways  to  fill  the  simulation  volume  with  different  media,  the  typical 
method  places  dielectric  media  wholly  within  a  cell  such  that  the  media  interfaces  are 
located  along  the  faces  of  the  Yee  grid.  In  this  way  the  continuity  of  the  tangential 
electric  fields  is  inherently  enforced.  The  construction  of  a  multi-dielectric  structure 
consists  of  assigning  a  permittivity  to  each  of  the  electric  field  component  positions 
within  the  computational  domain.  The  permittivities  are  based  on  the  physical  di¬ 
mensions  of  the  structure.  For  those  electric  field  components  that  lie  completely 
within  a  dielectric  medium  with  permittivity  Ci  the  permittivity  is  set  to  ei.  However 
for  those  electric  field  components  that  lie  on  the  planar  interface  between  two  media 
with  ei  and  €2,  an  average  of  the  two  dielectric’s  permittivities  is  used:  Finally, 

for  those  electric  field  components  that  lie  along  a  line  that  lies  at  the  intersection 
of  four  dielectric  media  with  ei,  £2,  ^3,  and  £4,  an  average  of  the  four  permittivities 
is  used;  Figure  2-4  is  a  diagram  that  shows  the  proper  assignment  of  the 

permittivity  values  for  a  FDTD  simulation  with  multiple  dielectric  type. 


2.6  Perfect  Conductors 


When  perfect  conductors  are  present  in  the  simulation  volume  the  the  total  electric 
field  must  be  zero  inside  the  conductor,  so  the  tangential  electric  fields  on  the  surface 
must  be  zero  to  satisfy  the  boundary  conditions.  If  an  analytic  plane  wave  excitation 
is  used  then  we  set 

E^  =  -E^  (2.40) 

inside  the  conductor. 

If  a  total  field  formulation  is  used  {i.e.  E^  =  0),  then  the  tangential  fields  are 
explicitly  set  to  zero  on  the  surface  of  the  conductor. 
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Figure  2-4:  Assignment  of  permittivity  values  for  multiple  dielectric  media  in  FDTD 
simulations. 


50 


CHAPTER  2.  FINITE  DIFFERENCE  TIME  DOMAIN  METHOD 


2.7  Accuracy,  Numerical  Dispersion,  and  Stability 

Since  the  FDTD  method  is  based  on  central  differencing,  the  smaller  the  grid  size 
the  more  accurately  the  difference  equations  will  approximate  the  partial  derivatives. 
However,  decreasing  the  grid  size  will  increase  the  memory  required  and  the  run  time 
of  the  simulation  so,  the  question  is  how  accurate  do  we  want  to  be? 

2.7.1  Accuracy  and  Numerical  Dispersion 

The  measure  of  accuracy  used  to  determine  the  grid  size  is  derived  from  the  numerical 
dispersion  relation  of  the  FDTD  grid  [14].  In  free  space  a  plane  wave  satisfies  the 
following  dispersion  relation: 


^  a;^eo)Uo, 


(2.41) 


where 


1 

—  —■ 

^0 

The  phase  velocity  in  any  direction  is  a  constant. 


U) 


'^ph  —  ^  —  Cq. 


(2.42) 


The  difference  between  the  numerical  dispersion  and  the  analytical  dispersion  is  the 
measure  of  the  accuracy.  The  dispersion  relation  of  the  FDTD  grid  with  a  uniform 
grid  spacing  is  given  by 


.  9  uAt 
sm  — — 
A 


(2.43) 


Instead  of  being  constant  the  FDTD  dispersion  relation  is  a  function  of  the  grid 
spacing.  A,  the  angle  of  propagation,  and  the  ratio,  called  the  stability  factor. 
As  expected,  when  At  0  and  A  — >  0,  equation  (2.43)  approaches  (2.41)  showing 
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Figure  2-5;  Normalized  phase  velocity  versus  angle  of  propagation  for  2D  FDTD;  ^  : 
(solid  line),  ^  :  (dashed  line),  and  |  :  (large-dashed  line). 

again  that  with  a  small  enough  grid  arbitrarily  good  accuracy  can  be  achieved  with 
FDTD.  Figure  2-5  is  a  plot  of  the  normalized  phase  velocity  versus  the  angle  that  the 
k  vector  makes  with  the  x  axis,  for  a  2D  FDTD  algorithm.  Each  curve  represents 
a  different  grid  size  in  terms  of  the  wavelength,  A.  The  solid  line  represents  the 
dashed  line  represents  ^  and  the  large-dashed  line  represents  |.  Linear  dispersion 
equates  to  a  normalized  phase  velocity  of  1  and  any  deviation  from  1  represents  error. 
The  FDTD  guideline  is  to  keep  the  dispersion  error  to  less  that  one  percent,  so  the 
grid  size  should  satisfy 

A  <  4-  (2-44) 


2.7.2  Stability 

Once  the  grid  size  has  been  determined  to  provide  a  reasonable  amount  of  accuracy, 
numerical  stability  places  an  upper  bound  on  the  time  step.  Stability  can  be  guaran- 
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teed  if  for  all  possible  wave  numbers  the  frequency,  u,  in  (2.41)  remains  real  valued. 
Thus,  the  stability  condition  for  a  non-uniform  cartesian  grid  is 


cAt  < 


1 


(AV)2 


+ 


~1~’ 

(A2)2 


(2.45) 


and  for  a  uniform  grid  is 

At  <  ^  (2.46) 

where  c  is  the  speed  of  light.  The  constraint  from  (2.45),  is  known  as  the  Courant 
stability  condition  and  must  be  enforced  at  all  times  within  all  materials. 

Many  FDTD  codes  will  allow  the  user  to  specify  the  grid  spacing  then  auto¬ 
matically  sets  the  time  step  to  the  maximum  allowed  under  the  stability  condition 
of  equation  2.46;  however,  sometimes  the  user  may  desire  a  smaller  time  step.  To 
accommodate  a  smaller  time  step  we  write 


At  < 


A 


CFLc 


(2.47) 


where  CFL  is  a  real  number  the  can  be  used  to  adjust  the  time  step,  keeping  in  mind 
that  there  is  a  lower  limit  and  that  limit  depends  on  the  dimension  of  the  simulation. 


2.8  Implementation  of  Sources  and  Excitations 


There  are  many  ways  to  excite  an  FDTD  simulation.  Whether  the  excitation  is  an 
analytic  plane  wave,  current  source,  or  boundary  condition,  a  Gaussian  pulse  is  used 
for  the  basic  temporal  excitation.  The  continuous  Gaussian  pulse  excitation  is 


g{t)  = 


(2.48) 
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The  Gaussian  pulse  is  desirable  because  it’s  spectrum  is  also  Gaussian  as  shown  below. 

OO 

G{uj)  —  j  g{t)U‘^*dt 

—  OO 

=  (2.49) 

Thus  a  Gaussian  pulse  may  be  used  to  excite  several  frequencies  in  a  single  simulation. 
In  FDTD,  since  the  Gaussian  pulse  must  be  truncated  {i.e.  the  signal  cannot  extend 
to  ±oo),  care  must  be  given  to  prevent  unwanted  high  frequency  components  from 
entering  the  simulation  from  the  waveform’s  truncation.  This  is  accomplished  by 
linking  L  and  T.  Discretizing  (2.48),  t  —>■  nAt  and  to  (3 At,  where  At  is  the  time 
step  and  n  and  /?  are  integers.  The  unwanted  frequencies  can  essentially  be  eliminated 
by  choosing  T  large  enough  to  damp  out  the  effects  of  the  truncation.  Typically, 

T=  to/4,  (2.50) 

so  that  the  attenuation  of  the  step  is  or  140  dB.  So,  the  pulse  width  of  the  signal, 
determined  by  to  =  PAt,  must  be  chosen  to  ensure  the  desired  bandwidth  and  T  is 
automatically  set  to  suppress  the  undesired  high  frequencies. 

The  pulse  width  of  the  excitation  is  approximately  2pAt  and  therefore  the  spec¬ 
tral  content  is  2^-  When  needed,  a  modulated  Gaussian  pulse  is  used  to  shift  the 
spectrum  to  a  desired  frequency  as  shown  below. 

g{t)  =  cos{u}ot)e~^^~*'°^^  (2-51) 

There  are  three  ways  to  insert  the  temporal  excitation  into  the  simulation.  The 
first  way  is  to  use  the  analytic  plane  wave  excitation,  E'‘,  in  equations  (2.35)-(2.37). 
The  second  way  is  to  add  electric  current  density  components  to  the  discretized 
Ampere’s  Law.  For  a  current  source  in  the  z  direction  the  cross  sectional  area  is 
assumed  to  be  Ax  Ay  and  a  length  of  Az  [20].  In  this  way  if  a  finer  grid  is  used 
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Figure  2-6:  2D  Yee  grid:  dashed  line  represents  the  edges  of  the  computational 
domain. 

for  the  same  physical  dimensions  then  the  current  density  component  must  be  scaled 
down  in  order  to  excite  with  the  same  current  density.  The  last  way  to  excite  the 
problem  is  with  a  boundary  condition  i.e.  define  an  electric  field  component  as  a 
function  of  time  at  some  location  within  the  simulation  volume. 


2.9  Treatment  of  Computational  Boundaries 

Figure  2-6  represents  a  two  dimensional  Yee  grid  for  a  TE  polarized  wave  with  fields 
Ey,Hx  and  H^.  The  dashed  lines  in  Figure  2-6  represent  the  truncation  of  the  com¬ 
putational  boundary  along  two  edges  to  the  domain’s  grid  lattice.  Since  the  magnetic 
field  quantities  above  and  to  the  right  of  the  dashed  lines  do  not  exist,  the  electric 
field  components  that  lie  along  the  boundaries  cannot  be  updated  with  the  standard 
difference  equations. 

When  the  edges  are  perfect  conductors,  the  tangential  electric  fields  can  be  set 
to  zero  and  therefore  the  boundaries  pose  no  difficulty.  However,  these  boundary 
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conditions  apply  to  a  narrow  range  of  electromagnetic  problems. 

If  left  alone,  the  edges  of  the  computational  domain  would  act  as  an  open  circuit 
and  reflect  the  outgoing  waves  that  impinge  upon  the  boundary.  With  a  large  enough 
domain,  measurements  can  be  taken  before  they  are  corrupted  by  the  non-physical 
reflections  from  the  boundary.  Unfortunately  this  technique,  called  time  gating,  may 
require  too  much  computer  memory  to  be  viable.  Many  approaches  have  been  used 
to  terminate  the  FDTD  computational  domain  [29]- [45]. 


2.9.1  One  Way  Wave  Equation 


The  most  common  treatment  of  the  boundaries  is  based  on  the  one  way  wave  equation 
made  popular  by  Mur  [29].  The  concept  is  simple:  construct  equations  that  allow 
only  outgoing  waves.  In  free  space  the  solutions  of  the  electric  field  will  satisfy  the 
wave  equation: 

+  +  =  (2-52) 


\dx‘‘  dy  dz  cq  dt 
Consider  a  boundary  at  a:  =  0  and  a  computational  domain  that  lies  to  the  right  of 
the  boundary,  a  plane  wave  traveling  in  free  space  is  given  by 


W{f,t)  = 


(2.53) 


where  W  represents  a  field  component  and 


k,  =  ±  (kl  -  kl  -  kl)  ' 


(2,54) 


with 


Re  S^(kl  -  kl  -  kiy^  >  0 


(2.55) 


1,2 

-  2  • 
Cq 


and 


(2.56) 
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Allowing  only  outgoing  waves,  then 


kx  =  -  {kl  ~kl~  kfj 


(2.57) 


The  next  step  is  to  convert  (2.57)  to  the  time  domain 


(2.58) 


where  Vy  and  are  the  phase  velocities  in  the  y  and  2  directions  defined  by 


Vy  = 


ui 


and 


= 


UJ 


The  essence  of  the  Mur  boundary  conditions  lies  in  the  approximation  to 


(l  -  {coVyf  -  {cov^f^  ■ 


(2.59) 


If  ((coUy)^  —  (cqUz)^)  1,  then  we  can  expand  (2.59)  as  follows; 

(l  -  {coVy)^  -  {cov^fy 
=  \-]^{{coVyf  F{coV^f) 

+  0{{{CoVyf  +  [cov^f)) 


The  first  order  Mur  uses  the  first  term  in  the  expansion.  This  is  equivalent  to 
assuming  normal  incidence  at  the  boundary.  The  partial  differential  equation  is 


A  _  il)  w 

.dx  Co  dt  I 


=  0 


(2.60) 


x—0 
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The  first  order  Mur  difference  equation  derived  from  (2.60)  is 


+  -  vr(0,3,k)].  (2.61) 


All  six  first  order  Mur  absorbing  boundary  conditions  can  be  found  in  Appendix  A. 
The  second  order  Mur  improves  the  approximation  by  using  the  first  two  terms 
in  the  expansion  so  that  the  partial  differential  equation  is 


1  d  d 
Co  dx  dt 


1 

cl  dt^  2 


(2.62) 


the  difference  equation  becomes 


[W^{Q,j,k)  +  W^{l,j,  k)] 


+ 


2Arr 


Co^i  "f"  Ax 


+l^”(l,i,  k  +  l)-  2W^{l,j,  k)  +  W”(l,  j,  k  -  1)] 

+W'"(1,3  +  1 1)  -  21V"(1,  j  +  1,  (:)  +  W“(l,  J  -  i.  fc)).  (2.63) 


All  twelve  second  order  Mur  absorbing  boundary  conditions  can  be  found  in  Ap¬ 
pendix  A. 


Coding  Considerations 

Since  the  Mur  ABC  equations  require  values  older  than  the  most  recent  time  step, 
the  Mur  equations  are  not  local  in  time  and  additional  memory  locations  are  neces- 
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Figure  2-7:  2D  Yee  grid  with  an  absorbing  boundary  layer  used  to  truncate  the 
computational  domain. 


sary  for  each  tangential  electric  field  on  the  boundary.  In  the  first  order  case,  the 
boundary  value  plus  its  nearest  neighbor  in  the  normal  direction  must  be  saved  at 
each  time  step  for  a  total  of  two  extra  memory  locations  per  tangential  field  compo¬ 
nent.  In  the  second  order  case,  two  past  time  steps  must  be  saved  for  a  total  of  four 
additional  memory  locations  per  tangential  field  component.  Although  higher  order 
equations  can  be  constructed,  the  additional  accuracy  does  not  usually  justify  the 
added  complexity  and  memory  required. 

For  efficient  implementation,  the  coefficients  found  in  the  Mm:  equations  are  cal¬ 
culated  once  at  the  beginning  of  the  program  and  stored  for  use  in  the  difference 
equations. 

2.9.2  Absorbing  Boundary  Layer 

In  addition  to  ABCs  based  on  one  way  wave  equation,  an  absorbing  layer  can  be  used 
to  absorb  the  waves.  A  layer  is  constructed  within  the  computational  domain  but 
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outside  the  simulation  volume.  The  idea  is  to  construct  a  layer  with  a  loss  mechanism 
that  will  allow  the  waves  to  propagate  into  the  absorbing  layer  with  little  reflection 
and  damp  out  the  waves  so  that  they  will  not  reflect  back  into  the  simulation  volume. 

Holland  [8]  showed  that  an  absorber  can  be  built  to  match  an  isotropic  lossless 
media  with  another  media  with  both  electric  and  magnetic  conductivities  as  long  as 
the  absorbing  media  has  the  same  permittivity  and  permeability  as  the  media  to  be 
matched  and 


a 

€ 


(2.64) 


With  (2.64)  set,  the  impedance  of  the  absorbing  layer  matches  the  impedance  of  the 
medium  inside  the  simulation  volume;  however,  this  type  of  ABC  is  only  good  at 
normal  incidence. 


2.10  Berenger’s  Perfectly  Matched  Layer 


In  1994,  Jean-Pierre  Berenger  introduced  a  novel  approach  [46]  to  the  absorbing 
layer  problem,  and  the  technique  has  been  the  subject  of  much  research  [47]-  [74]. 
Berenger  essentially  built  an  artificial  electromagnetic  absorbing  layer  around  the 
FDTD  computational  domain  that  offers  very  small  reflections  and  dissipates  the 
electromagnetic  energy  as  the  waves  travel  through  the  layer. 

Berenger’s  goal  was  to  remove  the  angle  dependence  of  Holland’s  ABC  [8]  while 
maintaining  the  reflectionless  property.  Berenger  modified  Maxwell’s  equations  by 
splitting  each  field  component  into  two  subcomponents,  one  for  each  of  the  other  two 
directions  as  shown  below: 

Ex  =  Exy  +  Exzi  (2.65) 

Ey  =  Eyx  +  Eyz,  (2.66) 

Ez  =  Ezx  +  Ezy,  (2.67) 


Hx  =  Hxy  +  H, 


(2.68) 
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Hy  =  Hyx  +  Hyz,  (2.69) 

and 

H,  -  H,x  +  H,y.  (2.70) 


He  further  added  anisotropic  electric  and  magnetic  conductivities.  Although  Ber- 
enger  originally  introduced  his  absorbing  boundary  for  two  dimensional  problems,  it 
was  quickly  expanded  to  three  and  the  12  equations  are 


^ ~  dy  >  (2-71) 

y^^Hxz  T  z  Hxz  —  {^yx  'E  Eyz) ,  (2.72) 

y'^Hyx  +  Cr^  Hyx  =  {Ezx  +  ^zy)  )  (2.73) 

~  ~^z  ^xz)  7  (2.74) 

d  d 

Hzx  —  ~‘q^  i.^yx  4"  7  (2.75) 

d  d 

^^y  ~  ^  i^xy  +  E^z)  j  (2.76) 

d  Q 

H“  ^yExy  ^  i^Hzx  "t"  Ezy^  ,  (2.77) 

d  d 

^"^Exz  “f"  ^zExZ  —  {Eyx  +  Hyz)  5  (2.78) 

^^^yx  +  ^xEyx  “  {Ezx  +  Hzy) ,  (2.79) 

4~  ^zEyz  ~  ~Q^  {Exy  -f^xz)  j  (2.80) 

^~^Ezx  +  cFj.Ezx  ^  {Eyx  +  Eyz)  5  (2.81) 

^'^^zy  +  *^y^zy  =  i^xy  +  Hxz)  ■  (2.82) 


and 
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The  key  to  understanding  the  PML  and  why  it  works  so  well  is  the  fact  that  if 
o-e  =  cr®  =  =  cr®  and  =  a'^  =  =  cr”*,  using  equations  (2.65)-(2.70),  the  PML 

equations  (2.71)-(2.82)  reduce  to  Maxwell’s  equations  of  a  lossy  medium  (2.18)  and 
(2.19).  In  other  words,  Maxwell’s  equations  are  a  special  case  of  the  more  general 
PML  equations. 

For  a  TE  polarized  wave  traveling  in  two  dimensional  space,  the  twelve  PML 
equations  reduce  to  six: 

+  alE,.  =  -L  (H,,  +  iJ,,)  (2.83) 


e-Ey,  +  alEy,  =  —  {H,y  +  H,,)  (2.84) 

l,^H,y  +  a^H,y  =  0  (2.85) 

^  {Ey,  +  Ey,)  (2.86) 

{Ey.  +  Ey,)  (2.87) 


txj^H^y  +  a^H^y  =  0  (2.88) 

Since  the  y  direction  collapses  and  =  0,  a  plane  wave  solution  yields 

Ey,  =  (2.89) 

uectl 

Ey,  =  (2.90) 

ueal 
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(2.91) 

(2.92) 

(2.93) 

(2.94) 

(2.95) 

(2.96) 

(2.97) 

(2.98) 

(2.99) 


Berenger  and  Cai  et  al  [51]  offer  analytical  derivations  of  the  PML  conditions  for 
matching.  The  following  is  another  that  offers  insight  into  the  PML  mechanism. 
Consider  two  PML  media  with  non-zero  values  of  permittivity,  permeability  and 
conductivity;  one  PML  medium  occupies  Region  0  and  the  other  Region  t.  A  TE 
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Figure  2-8;  Geometry  for  PML  reflection  coefficient  analysis. 


polarized  incident  wave  with  k  =  k^x  —  k^z,  strikes  the  boundary  from  above.  The 
electric  fields  are  assumed  [137]  to  be  of  the  form 


Ei  =  y{Ey,o  +  Ey,o)e^^^^-^'‘^^ 


(2.100) 


Er  =  yR'^^  {Ey^o  +  Ey,o)  (2.101) 

and 

Et  =  yT^^  {Ey^o  +  Ey,o)  (2.102) 

Using  (2.91)  and  (2.92)  to  determine  the  magnetic  fields  and  imposing  the  boundary 
conditions  at  2:  =  0  the  reflection  coefficient  is 
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The  reflection  characteristics  can  be  found  in  the  term 


(2.104) 


From  the  PML  dispersion  relations  for  each  PML  medium  (2.104)  becomes 


a^at^a%af  -  kl^) 

\  OCuOitO^txOiZ  -  kli) 


(2.105) 


For  the  interface  to  be  reflectionless  at  all  angles,  the  k  dependence  must  be  removed. 
This  is  possible  if 

=  oclafk,  (2.106) 

which  is  true  if  et  =  e  ,  fit  =  fJ.  and  angle  dependence  is  removed 

leaving 

ln/mn/<r 

(2.107) 


pIF 


a%aZ  ' 


Finally  if  Poi  =  !>  interface  will  be  reflectionless  which  is  satisfled  by  must 

equal  ala^.  In  terms  of  the  media  parameters,  to  match  {i.e.  reflectionless  at  all 
angles)  the  PML  medium  in  Region  0  with  (e,  /i,  <7®,  a^,  cr®,  cr^)  with  the  PML  media 
in  Region  t,  Region  fs  parameters  must  be 


Ct  —  e. 

(2.108) 

Q  =  e. 

(2.109) 

II 

(2.110) 

II 

b 

(2.111) 

b 

II 

b 

(2.112) 
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_m  O'  e 

=  -^z- 
€ 


(2.113) 


The  last  condition  (2.113)  is  imposed  on  Region  0.  Thus,  the  electric  and  magnetic 
conductivities  in  Region  0  can  not  be  arbitrary.  However,  a  lossless  medium  with  e, 
fx,  does  satisfy  (2.113),  and  the  corresponding  matching  parameters  in  Region  t  are 


€t  =  e, 


(2.114) 


Ot  =  0, 


=  0, 


_  r\ 

^tx  b. 


(2.115) 

(2.116) 
(2.117) 


_m  _  0  e 

^tz  ~  ^  '^tz- 


(2.118) 


For  infinite  half-spaces,  the  selection  of  is  arbitrary.  In  FDTD  simulations,  how¬ 
ever,  the  selection  of  is  not  arbitrary  and  will  be  discussed  in  a  future  section. 


2.10.2  PML  Propagation  Characteristics 

Again  consider  the  transmission  from  a  lossless  dielectric  in  Region  0  into  its  matched 
PML  in  Region  t.  From  the  dispersion  relation  in  Region  t, 


-LLL-  +  =  kf, 


(2.119) 


and  the  incident  k  vector  x-component,  the  dispersion  relation  in  Region  t,  sim¬ 


plifies  to 


7.2  ,  '^tz  _  7  2 


(2.120) 
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Therefore,  the  z  component  of  the  k  vectors  are  related  by  ktz  =  Otz^iz  or 

ktz  =  kiz  +  i—kiz.  (2.121) 

cue 

When  the  wave  enters  the  PML,  it  propagates  in  the  same  direction  as  the  incident 
wave  and  is  attenuated  in  the  direction  normal  to  the  interface.  In  terms  of  the 
incident  angle,  (2.121)  becomes 


ktz  ^  kiz  +  cos  9.  (2.122) 

The  imaginary  part  of  the  ktz  is  a  function  of  the  incident  angle,  providing  the  most 
damping  in  the  normal  direction  and  no  damping  in  the  tangential  direction.  If 
the  incident  wave  is  evanescent  (  kiz  purely  imaginary),  then  the  PML  provides  no 
additional  damping. 

2.10.3  PML  and  FDTD  Simulations 

So  far,  only  infinite  half-spaces  have  been  considered.  An  infinite  half-space  is  not  ap¬ 
propriate  for  FDTD  simulations  since  the  computational  domain  must  be  truncated 
at  some  point.  Berenger  proposed  truncating  the  PML  with  a  Perfect  Electric  Con¬ 
ductor  (PEC)  i.e.  setting  all  tangential  fields  to  zero  at  the  outer  boundary.  Thus 
a  wave  will  travel  through  the  PML,  reflect  off  of  the  PEC,  travel  through  the  PML 
again,  and  re-enter  the  simulation  volume.  However,  the  wave  will  be  attenuated 
by  the  loss  mechanism  (2.121)  in  the  PML.  The  user  will  select  a  PML  thickness,  a 
PML  conductivity,  and  a  conductivity  profile  to  achieve  the  desired  low  reflection  co¬ 
efficient.  Even  though  analytic  analysis  shows  that  a  constant  conductivity  profile  is 
sufficient  for  a  zero  reflection  at  the  boundary,  numerical  reflections  will  occur  at  any 
interface  that  contains  conductivity  discontinuities [60].  To  reduce  these  numerical  re¬ 
flections,  Berenger  proposed  grading  the  conductivity  from  a  low  value  of  zero  at  the 
boundary  to  some  maximum  value  at  the  PEC  wall.  The  most  common  conductivity 


2. 1 0.  BERENGER  ’S  PERFECTLY  MATCHED  LAYER 


67 


profile  is  given  by 

(9”  (2  123) 

where  p  is  the  normal  coordinate  variable,  d  is  the  thickness  of  the  PML,  and  n  is 
the  conductivity  grading  factor  so  that  n  =  0  is  a  constant  profile,  n  =  1  is  a  linear 
profile,  n  =  2  is  a  parabolic  profile,  etc.  If  the  PML  starts  at  p  =  0  and  ends  at  p  =  d 
then  the  expression  for  the  reflection  coefficient  is  given  by 

d 

-2  f  a^(p)7j  cos  0  dp 

R(0)  =  e  0  ,  (2.124) 


and  after  the  integration 


-2tj  da-' 


R(9)  =  e — ^ 


P  ’^9:7.  cos  e 


(2.125) 


So  for  practical  FDTD  applications  the  reflection  coefficient  is  still  a  function  of  the 
incident  angle;  however,  the  low  grazing  angle  waves  will  be  absorbed  by  the  PML 
that  terminates  the  simulation  volume  at  the  edges  perpendicular  to  the  normal  PML 
layer.  Also,  it  appears  that  the  reflection  coefficient  can  be  made  arbitrarily  small  with 
a  large  enough  conductivity,  but  as  mentioned,  numerical  reflections  occur  with  large 
changes  in  conductivity  between  FDTD  cells  and  prevents  arbitrarily  small  reflection 
coefficients. 

Berenger  proposed  that  one  way  to  select  the  conductivity  is  to  set  the  reflection 
coefficient  at  normal  incidence  to  some  desired  value.  Setting  6  =  0  and  solving  for 
(Tpmaxi  conductivity  is  determined  by 


'  p  max 


{n  +  1)  In  i?  (0) 
2dri 


(2.126) 


2.10.4  PML  and  Multilayer  Dielectrics 

When  applying  Berenger’s  PML  to  multilayer  dielectrics,  as  in  the  case  of  microstrip 
and  dielectric  waveguides,  special  care  must  be  taken  when  determining  the  conduc- 
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Figure  2-9:  Configuration  for  multi-dielectric  layer  PML. 


tivities  for  each  layer.  For  a  single  layer  {i.e.  no  layer)  the  conductivity  is  calculated 
based  on  the  desired  reflection  coefficient  at  normal  incidence. 


Now  consider  two  lossless  dielectrics  with  permittivities,  ei  and  62,  and  permeabil¬ 
ity,  //(,,  in  Region  0  with  their  interface  at  x  =  0.  Layer  1  is  located  at  x  <  0  and  layer 
2  is  located  at  x  >  0.  The  layers  are  terminated  with  a  PML  in  Region  t  with  two 
conductivities,  an  and  <742.  See  Figure  2-9.  Given  the  same  PML  thickness,  profile, 
and  desired  reflection  coefficient,  if  computed  independently,  the  conductivities  are 


a 


e 

Itz  max 


{n  -f  1)  InR(O) 
2drji 


(2.127) 


and 


<7 


e  _ 

2tzmax 


(n+l)lnR(0) 

2dri2 


(2.128) 
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In  this  way,  the  conductivity  ratio  between  the  two  PML  terminating  layers  is 


^  2tz  max 
^Itz  max 


(2.129) 


However,  if  the  ratio  of  the  reflection  coefficients  is  examined;  given  by 


Ron  {0) 

R(St2  (^) 


-2dy/fio{ 


1  mag 


_ 211 _ 


=  e 


n+1 


(2.130) 


Substitution  of  (2.129)  into  (2.130)  leads  to 

Ron  (^)  — COS02) 

Rot2{0)~ 


(2.131) 


The  reflection  coefficients  will  be  equal  only  when  9i  =  62  which  is  true  only  at  grazing 
incidence  or  when  there  is  no  boundary  at  x  =  0.  Both  cases  are  uninteresting.  On  the 
other  hand,  if  phase  matching  in  Region  0  is  applied  x  =  0  i.e.  A;icos0i  =  k2Cos92, 
substitution  leads  to 


-Rpti  {9) 
Rot2  (9) 


—  2dy/ixo{ 


”1  max  ”2  mnx2^ 

Vn _ "2 


)  COS  9i 


—  e 


n+1 


(2.132) 


In  order  for  the  reflection  coefficients  along  2;  =  0  to  be  continuous,  i.e. 


Ron  (^)  _  , 

Rot2  (9) 

the  conductivities  between  adjacent  dielectric  layers  must  satisfy 


(2.133) 


e  _ 

^2max  ~~  ^ 


1  max' 


(2.134) 


Further  insight  into  this  condition  can  be  seen  from  examining  the  reflection  coeffi¬ 
cients  in  the  direction  tangential  to  the  dielectric-PML  interface,  at  the  layer  1/layer  2 
interface  at  a;  =  0.  The  reflection  coefficient  at  the  dielectric  interface,  Rji^iz  >  0), 
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Figure  2-10:  Configuration  for  multi-dielectric  reflection  coefficient  analysis. 


should  equal  the  reflection  coefficient  at  the  matching  PML  layer  interface,  Rj-fiz  < 
0),  for  the  PML  to  truly  simulate  open  space.  Prom  the  original  analysis  of  the  PML, 
the  reflection  coefficient  between  any  two  PML  media  is 

Rif  =  ^  (2.135) 


where 


pli  = 


\  {a\^aZkl  -  kf) ' 


(2.136) 


After  the  setting  the  material  parameters  and  PML  matching  conditions  between 
Regions  0  and  t  ,  (2.136)  becomes 


U2  _  h,2 
u2  _  h.2 


Pniz>0)  = 


(2.137) 
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and 


pjiiz<0)  = 


C^ltzO^Uz  -  k%) 

^  {o^itzO^uz^it  -  kit) ' 


(2.138) 


Phase  matching  at  the  x  =  0  interface  says  ku  =  ^22  =  ^2  aiid  kuz  =  k2tz  =  hz-  The 
PML  matching  conditions  along  with  the  PML  dispersion  relation  imply  ku  =  ki, 
k2t  =  k2,  a|t2  =  ^2tz^  kuz  =  altz^z  and  k2tz  =  Substitution  of 

these  relationships  into  (2.138)  reveal 


pff  {z<0)  = 


\  kl-kl 


(2.139) 


Equations  (2.138)  and  (2.138)  are  equivalent  if  an  additional  PML  condition  is  met: 


a 


e 

Itz 


2tz- 


(2.140) 


This  additional  condition  is  equivalent  to  (2.64).  Ensuring  that  the  reflection  coef¬ 
ficients  be  continuous  at  the  intersection  of  the  four  different  media  leads  directly 
to 


*^2  max  ^ 


e2_e 

1  max’ 

^1 


(2.141) 


Enforcing  this  ratio  is  equivalent  to  saying  the  loss  tangents  within  each  PML  layer 
must  be  equal.  Bahr  et  al  in  [50]  arrived  at  the  same  conclusion  by  demanding  that 
there  be  equal  decay  in  the  normal  direction  in  each  of  the  layers  with  phase  matching 
in  mind. 


Failure  to  impose  this  condition  can  lead  to  instabilities  at  the  PML/PML  inter¬ 
face  due  to  the  generation  of  other  waves  necessary  to  match  the  boundary  conditions. 
In  the  next  chapter,  the  generation  of  mode  templates  for  dielectric  waveguides  is  dis¬ 
cussed.  When  generating  these  templates  many  time  steps  are  required  to  capture 
the  necessary  frequency  resolution.  When  a  large  permittivity  contrast  between  ex¬ 
ists  between  two  layers  in  a  multi-layer  dielectric  structure,  the  independent  PML 
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Figure  2-11:  Test  setup  for  PML  instability  demonstration.  Field  measured  at 
point  A. 


matching  condition  (2.129)  produces  instabilities. 

For  example,  the  fundamental  mode  spatial  distribution  is  calculated  from  a  2D 
FDTD  simulation  of  a  dielectric  rib  built  with  Gallium  Arsenide  with  €2  =  ll.Sco 
on  top  of  a  substrate  €3  =  lO.Oeo  surrounded  with  free  space  cladding,  cj  =  cq.  See 
Figure  2-11.  The  grid  size  is  A  =  .002  m  and  the  time  step  is  set  to  the  3D  Courant 
limit.  The  simulation  volume  is  57  x  41. 

Figure  2-12  represents  E^  versus  time  step  at  Point  A  in  Figure  2-11  within  the 
FDTD  simulation  volume  with  two  different  PML  matching  conditions.  Both  sim¬ 
ulations  use  an  eight  layer  PML  with  n  =  2  and  i?(0)  =  10“®  to  terminate  the 
computational  domain.  The  conductivities  of  each  layer  were  calculated  using  the 
dependent  matching  conditions  of  equation  2.134  and  independent  matching  condi¬ 
tions  of  equation  2.129.  Table  2.1  provides  the  conductivity  values.  The  two  graphs 
clearly  show  the  effects  of  using  the  independent  matching  condition.  The  indepen¬ 
dent  construction  of  the  PML  conductivities  produces  solutions  that  grow  and  take 
over  the  simulation  after  3000  time  steps. 
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Figure  2-12:  Field  Ez  measurements  at  PML/PML  interface  (Point  A  in  Figure  2-11) 
with  dependent  (top)  and  independent  (lower)  PML  matching  conditions. 
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Layer 

Permittivity 

e 

^max 

Dependent 

Independent 

Cladding 

l.Oeo 

11.45 

11.45 

Film 

ll.Seo 

135.11 

39.33 

Substrate 

lO.Oeo 

114.50 

36.21 

Table  2.1:  Maximum  conductivities  to  match  PML  to  rib  structure. 

The  simulations  were  constructed  to  excite  the  fundamental  mode  of  the  rib  waveg¬ 
uide  by  fixing  the  propagation  constant  at  600  and  exciting  the  problem  with  a 
wide  {P  —  256)  Gaussian  pulse  modulated  at  the  corresponding  mode  temporal  fre¬ 
quency  of  /o  =  8.1  GHz.  Figure  2-13  represents  the  temporal  frequency  content  inside 
the  rib  waveguide  due  to  the  modulated  Gaussian  pulse  with  with  a  fixed  propagation 
constant  with  the  dependent  (top)  and  independent  (lower)  PML  matching  conditions. 
The  lower  spectrum  shows  the  effects  of  the  independent  PML  matching  conditions. 
The  singularities  at  the  PML/PML  interfaces  act  as  sources  and  place  energy  in 
frequencies  other  than  the  fundamental  mode.  These  sources  make  it  difficult  to  de¬ 
termine  the  fundamental  mode  spatial  arrangement  as  shown  in  Figure  2-14.  The 
lower  plot  has  unwanted  spikes  in  the  mode’s  spatial  distribution  at  the  PML/PML 
interfaces  as  seen  on  the  left  and  right  sides  of  the  computational  domain. 

2.10.5  PML  Coefficients  in  FDTD 

The  PML  difference  equations  use  an  explicit  exponential  form  as  proposed  by  Ber- 
enger.  For  example,  the  E^y  equation  used  in  the  PML  region  of  the  computational 
domain  is  given  by 

+  =  C(’li  +  i,j,k)B;^{i  +  i,j,k)  (2.142) 

-  +  +  kj  +  k,k)  -  +  kj  -  k,k) 
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Figure  2-13:  Fourier  transform  of  field  inside  rib  using  dependent  (top)  and  indepen¬ 
dent  (bottom)  PML  conductivity  matching  condition. 
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Figure  2-14:  Fundamental  mode  spatial  field  distribution  using  dependent  (top)  and 
independent  (bottom)  PML  conductivity  matching  condition. 


2.11.  CONSIDERATIONS  AND  OVERVIEW  OF  FDTD  ALGORITHM 
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where  the  two  constants,  and  are  given  by 


(2.143) 


and 


/j  _L  i  ^  lr\  —  1 - 1 - 


cr{i+ij^)=  ,.^1  .... 


(2.144) 

n  ic  \  /\n 

'y\^  '  2’ 

The  remainder  of  the  PML  FDTD  difference  equations  along  with  a  very  detailed 
description  of  the  PML  coefficients  is  found  in  Appendix  B 


2.11  Considerations  and  Overview  of  FDTD  Algo¬ 
rithm 

The  FDTD  method  is  a  very  versatile  method  for  solving  electromagnetic  problems. 
The  discretization  of  Maxwell’s  equations  has  a  built  in  error  that  is  proportional 
to  the  square  of  the  grid  size  (A^).  In  addition,  the  difference  equations  introduce 
non-physical  numerical  dispersion.  The  dispersion  error  can  be  kept  to  less  than 
one  percent  if  the  grid  size  is  less  than  one  tenth  of  a  wavelength.  When  there  is 
more  than  one  material  in  the  simulation  volume,  it  is  the  most  dense  material  that 
determines  the  grid  size.  This  means  that  the  less  dense  material  is  essentially  over 
sampled.  If  two  grid  sizes  are  used,  then  numerical  reflections  may  result  [14].  Also, 
the  FDTD  method  must  model  curved  surfaces  with  a  staircase.  When  the  other 
errors  are  kept  small,  the  grid  size  is  typically  sufficiently  fine  that  the  staircase  is  a 
good  approximation  to  a  curved  surface.  Finally,  there  will  always  be  some  reflections 
from  the  boundaries  of  the  computational  domain  no  matter  which  ABC  is  used  and 
must  be  considered  when  constructing  the  problem. 

The  basic  steps  to  the  FDTD  algorithm  are  listed  below. 


1.  Build  geometry. 
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2.  Initialize  variables  and  calculate  coefficients. 

3.  Calculate  electric  fields. 

4.  Impose  source  condition. 

5.  Advance  time  by  half  step. 

6.  Calculate  magnetic  fields. 

7.  Advance  time  by  half  step. 

8.  Save  desired  values. 

9.  Repeat  time  marching  until  steady  state  is  achieved. 

10.  Perform  post-processing  on  saved  data. 

2.12  Summary 

In  this  chapter  the  Finite  Difference  Time  Domain  method  for  solving  Maxwell’s  equa¬ 
tions  has  been  presented.  The  method  is  simply  a  finite  difference  approximation  to 
the  partial  derivatives  found  in  Maxwell’s  curl  equations.  Electromagnetic  wave  in¬ 
teractions  with  fairly  complicated  structures  can  be  modeled  and  multiple  frequencies 
can  be  tested  by  controlling  the  pulse  width  of  the  excitation. 

Difference  equations  suitable  for  coding  were  presented  for  free  space,  lossy  dielec¬ 
tric  media,  Mur  absorbing  boundaries,  and  the  Perfectly  Matched  Layer. 

The  artificial  PML  conductivities  used  to  absorb  the  outgoing  waves  in  a  multi¬ 
ple  dielectric  layer  simulation  must  not  be  assigned  independently  from  Berenger’s 
proposed  method  of  normal  incidence;  instead,  once  one  of  the  layer’s  conductivities 
is  calculated  the  others  must  be  constructed  so  that  the  loss  tangents  in  each  layer 
are  equal.  With  this  done,  there  will  be  no  reflection  coefficient  singularities  in  both 
the  normal  and  tangential  directions  and  the  PML  will  simulate  open  space  to  the 
maximum  extent  possible  under  the  discretization  scheme. 


Chapter  3 


FDTD  Analysis  of  Dielectric  Rib 
Waveguide  Discontinuities 

3.1  Introduction 

Dielectric  waveguides  are  commonly  used  as  interconnects  in  millimeter-wave  and 
sub-millimeter  wave  integrated  circuit  technologies  [75]- [100].  These  interconnects 
contain  discontinuities  such  as  bends  that  can  introduce  high  loss  through  the  waveg¬ 
uide.  Since  the  propagating  modes  of  this  type  of  structure  are  characterized  by  very 
complex  field  distributions,  with  no  exact  analytic  solutions,  a  numerical  approach  is 
appropriate  to  investigate  these  structures.  In  this  chapter,  a  full  three  dimensional 
FDTD  implementation  is  used  to  analyze  the  effects  of  waveguide  discontinuities  on 
the  fundamental  mode’s  propagation  through  the  waveguide. 

When  using  FDTD  to  study  any  waveguide,  depending  on  the  excitation  used,  it 
may  take  time  for  the  desired  mode  {e.g.  the  fundamental)  to  develop  as  the  evanes¬ 
cent  modes  decay  away  and  radiating  modes  leave  the  simulation  volume.  In  order  to 
launch  the  fundamental  mode  more  efficiently,  in  this  work  the  mode’s  spatial  distri¬ 
bution  is  determined  up  front.  A  two  dimensional  FDTD  method  is  used  to  construct 
a  mode  template  [131,  133].  The  three  dimensional  FDTD  code  is  collapsed  in  the 
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direction  of  propagation  by  replacing  the  spatial  difference  equation  with  an  equa¬ 
tion  based  an  assumed  propagation  constant  and  the  corresponding  phase  difference 
between  adjacent  spatial  grid  points.  A  two  dimensional  simulation  of  the  waveguide 
cross  section  produces  the  corresponding  temporal  frequencies  of  the  propagating 
modes.  The  lowest  frequency  is  the  fundamental  mode  and  the  spatial  distribution  is 
found  by  performing  a  Fourier  transform  at  each  space  point  on  the  two  dimensional 
grid  at  the  fundamental  mode  frequency.  With  this  source  condition,  the  compu¬ 
tational  domain  can  be  reduced  from  the  commonly  used  spatial  Gaussian  source 
condition  by  allowing  shorter  distances  between  excitation  and  the  discontinuity. 

The  FDTD  method  has  been  used  to  study  metallic  waveguides  and  metallic  pla¬ 
nar  structures  [101]-[108].  Unlike  metallic  waveguides  the  dielectric  rib  waveguide 
is  an  open  structure  and  analyzing  them  with  the  Finite  Difference  Time  Domain 
method  requires  absorbing  boundary  conditions  that  simulate  open  space.  The  per¬ 
fectly  matched  layer  (PML)  described  and  analyzed  in  Chapter  2  is  used  to  terminate 
the  simulation  volume  unless  otherwise  stated. 

With  mode  template  excitations  and  PML  absorbing  boundary  conditions,  the 
dielectric  rib  waveguide  is  studied.  Different  waveguide  bends  are  examined  along 
with  the  effects  of  the  rib  geometry. 

3.2  Implementation  of  Dielectric  Rib  Waveguide 
in  FDTD 

Various  dielectric  rib  waveguides  are  constructed  and  analyzed  in  this  chapter.  The 
input  parameters  are  the  rib  width  (rw),  rib  height  (r/i),  film  height  (fh)  and  the 
center  reference  point  {xc,zc)  as  shown  in  Figure  3-1.  The  simulation  model  of  a 
waveguide  is  constructed  by  carefully  assigning  a  permittivity  at  each  electric  field 
location  as  described  in  Section  2.5.  The  initial  propagation  is  in  the  +y  direction. 

The  Berenger  perfectly  matched  layer  described  in  Section  2.10  is  used  to  termi- 
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Figure  3-1:  Schematic  of  dielectric  rib  waveguide  cross  section.  Parameters  are  rw  : 
rib  width,  rh  :  rib  height,  fh  :  film  height,  and  {xc,  zc)  :  center  reference  point. 

nate  the  FDTD  computational  domain  for  all  three  dimensional  simulations  unless 
otherwise  stated.  The  PML  consists  of  8  discretized  layers  using  a  quadratic  con¬ 
ductivity  profile  {i.e.  n  =  2)  with  the  maximum  conductivity  derived  by  setting 
R{0)  =  10”^  with  free  space  permittivity,  eo-  The  conductivities  for  the  other  dielec¬ 
tric  layers  are  calculated  with  equation  (2.134). 

To  analyze  waveguide  discontinuities,  the  S  parameters  are  calculated  from  the 
FDTD  simulation  measurements  of  the  time-domain  waveforms.  Multiple  frequen¬ 
cies  are  analyzed  within  a  single  FDTD  simulation.  The  methods  are  described  in 
Section  3.7.1. 


3.3  Approximate  Analytic  Solution 

The  study  of  dielectric  waveguide  modes  and  their  propagation  characteristics  have 
been  studied  since  the  late  sixties,  except  for  the  dielectric  slab  and  circular  dielec¬ 
tric  waveguides,  exact  analytical  solutions  do  not  exist  and  numerical  techniques  are 
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In  1969,  Marcatili  introduced  an  approximate  analytical  solution  for  the  guided 
modes  of  a  rectangular  dielectric  waveguide  that  is  still  used  today.  Figure  3-2 
represents  the  cross-section  of  a  rectangular  dielectric  waveguide  with  permittivity, 
ei,  surrounded  by  four  different  dielectric  media,  62,3,4,5  with  dispersion  relations: 
kf  =  uj^nei,  Z  =  1  —  5.  As  summarized  in  [142],  Marcatili  derived  an  approximate  so¬ 
lution  by  applying  the  boundary  conditions  to  the  edges  of  the  guide  only.  He  assumed 
sinusoidal  variations  along  the  x  and  2:  directions  inside  the  guide  and  exponential 
decay  outside  of  the  guide. 

For  the  ^  modes  the  governing  equations  are 

k^a  =  SIT  -  tan"^(A;^^3)  -  tsin~^{kxC5),  (3.1) 

and 

k^b  =  qir  —  tan“^  ^  ^2  2) 
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where 


6,5  — 


C2,4  —  — 


Xi  =  n[ki-kf)U  =  2-5.  (3.5) 

Phase  matching  at  2;  =  0  and  z  =  b  implies  =  kxi  =  kx2  =  and  similarly  at 
X  =  0  and  x  =  a  implies  kz  =  hi  =  k^s  =  k^s-  Propagation  is  in  the  y  direction  with 
a  propagation  constant  given  by 


I3^ky^  {kl-kl-  kfj " . 


For  the  dielectric  waveguide  surrounded  by  the  same  dielectric  such  as  free  space, 
62  =  63  =  64  =  €5  =  Co,  the  equations  simplify  to 

fcajU  =  STT  —  2tan“^(A:x(^),  (3.7) 


hh  =  q-K  —  2  tan  ^  ^ 


6,5  —  ^ 


C2,4  =  C  = 


(3.10) 


X2,3,4,5  —  X  — 


(3.11) 
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At  a  given  frequency,  oj,  x  can  be  found  with  equations  (3.7)  and  (3.11).  Substitution 
of  X  iiifc  (3.9)  and  (3.10)  give  ^  and  Graphical  techniques  or  Newton’s  method 
can  be  used  to  find  the  discrete  set  of  and  kz  using  (3.7)  and  (3.8)  that  correspond 
to  a  given  mode  number  and  uj.  However,  the  Marcatili  Method  is  only  good  when 
most  of  the  mode’s  energy  is  confined  to  the  core  dielectric  region  and  is  not  reliable 
near  cutoff  [76]. 

Dispersion  curves  for  dielectric  waveguides  are  often  displayed  in  terms  of  normal¬ 
ized  variables.  The  normalized  propagation  constant  is  defined  by 


and  the  normalized  waveguide  height  is 


(3.12) 


(3.13) 


where  A  =  ^  and  Aq  = 

In  order  to  validate  the  upcoming  FDTD  mode  construction  procedure,  a  square 
dielectric  waveguide  is  examined.  The  waveguide  has  a  permittivity  of  2.8eo  sur¬ 
rounded  by  a  cladding  of  free  space.  The  structure  measures  1  cm  on  a  side.  Using 
equations  (3.7)-(3.11),  the  dispersion  curves  were  constructed  for  the  F'f  i,  ^^^d 
£^^21  modes.  Figure  3-3  represents  the  normalized  curves  using  (3.12)  and  (3.13).  The 
results  of  the  fundamental  mode  will  be  used  to  compare  with  the  FDTD  method. 


3.4  Excitation 

The  excitation  consists  of  a  source  plane  of  distributed  z  directed  current  sources  at 
the  interface  between  the  simulation  volume  and  the  PML  that  forms  the  front  of  the 
computational  domain.  This  excitation  technique  has  advantages  over  the  an  electric 
field  boundary  condition;  a  major  one  being  that  the  PML  can  absorb  waves  reflected 
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First  Three  Modes 


Figure  3-3;  Normalized  dispersion  curves  of  first  three  modes,  (diamond)  £'|,i) 
(triangle)  and  Ef  2  (star),  as  derived  by  the  Marcatili  method. 


from  the  waveguide  discontinuities  [50] .  When  Mur  boundary  conditions  are  used  one 
of  the  faces  has  a  dual  role:  first  as  an  excitation  plane  and  then  as  an  absorbing 
boundary.  After  the  excitation  is  completely  launched,  the  source  is  essentially  turned 
off  and  the  absorbing  boundary  condition  is  turned  on  [106].  Using  this  method,  a 
larger  computational  domain  is  necessary  so  that  the  ABC  can  be  turned  on  after 
the  excitation  but  before  any  reflections  from  the  discontinuity  reach  the  source  edge. 
The  magnitude  of  the  individual  current  sources  are  derived  up  front  with  a  template 
calculation  described  in  Section  3.5.  The  temporal  excitation  is  that  of  a  modulated 
Gaussian  as  described  in  Section  2.8. 
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3.5  Improvement  to  the  FDTD  Method  with  a 
Better  Excitation 

Although  the  recent  strides  in  computer  technology  have  placed  more  memory  and 
computing  power  in  the  hands  of  engineers,  reducing  computational  domain  require¬ 
ments  is  always  desirable.  The  reduction  in  the  domain  leads  to  faster  simulations 
and  the  ability  to  simulate  larger  problems.  In  this  section  the  the  construction  of  a 
mode  template  is  described  and  tested  to  show  how  the  computational  domain  can 
be  reduced. 

3.5.1  Template  Construction 

When  performing  discontinuity  analysis,  it  is  desirable  to  launch  a  known  mode.  How¬ 
ever,  complicated  structures  like  the  rib  waveguide  have  modes  that  have  complicated 
spatial  distributions.  These  modes  can  not  be  expressed  analytically;  therefore,  a  suf¬ 
ficient  length  of  section  of  computational  domain  after  the  source  plane  and  before 
the  discontinuity  is  needed  to  allow  the  fundamental  mode  to  set  up. 

A  numerical  approach  has  been  used  with  closed  and  partially  closed  structures  to 
find  what  is  called  a  mode  template  [131, 133].  The  template  is  the  spatial  arrangement 
of  the  field  in  the  plane  perpendicular  to  the  direction  of  propagation.  In  this  section, 
the  numerical  technique  of  templates  is  extended  to  open  structures. 

To  construct  a  mode  template,  the  three  dimensional  FDTD  algorithm  is  con¬ 
verted  to  a  two  dimensional  algorithm  where  the  propagating  direction  is  replaced 
by  equations  that  incorporate  the  assumed  known  phase  difference  between  planes 
perpendicular  to  the  propagating  direction. 

The  first  step  is  to  assume  a  propagation  constant,  ky.  Then  the  y  dependence  in 
the  difference  equations  is  modified  by 


k)  -  F{i,j  -l,k)^  F{i,j,  k)(l  -  e-'*"**) 


(314) 
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and 

F(2j  +  l,k)-  F(i,j,  k)  ^  F{i,j,  k){e^^y^y  -  1).  (3.15) 

With  the  new  algorithm  and  a  propagation  constant,  a  two  dimensional  simulation  is 
run  and  the  time  domain  waveform  is  saved.  The  Fourier  transform  is  taken  of  the 
steady  state  portion  of  signal.  The  first  half  of  the  signal  is  discarded.  The  peaks  in 
the  Fourier  transform  correspond  the  mode  frequencies  associated  with  the  selected 
propagation  constant. 

The  rationale  behind  this  template  construction  can  be  seen  in  the  following  anal¬ 
ysis.  We  assume  that  the  steady  state  of  a  component  of  the  electric  field  can  be 
expressed  as  the  sum  of  propagating  modes: 


Ea(x,  y,  Z.  t)  =  yi  z)e'‘^''e  a  =  x,y,  z. 


(3.16) 


The  Fourier  transform  of  the  electric  field  is 


TOO 

Ea{x,y,z,uj)  =  /  Ea{x,y,z,t)F‘^Mt,  a  =  x,y,z. 

J  —  OO 


(3.17) 


Substitution  of  (3.16)  into  (3.17)  yields 


roo 

B„(i,!/,z.u,)  =  /  y0„(x,z)e-‘«''e-""‘e'“‘<ft 

J-oo  n 

/OO 

-°o 

=  4>n{x,  z)e"^y^2'K6{uJ  -  U>n) 

n 

a  =  x,y,z. 


(3.18) 


The  delta  functions  in  (3.18)  implies  that  peaks  in  the  spectrum  will  correspond  to 
the  mode  frequencies. 

Once  the  mode  frequency  is  determined,  the  Fourier  transform  of  each  {x,  z)  lo¬ 
cation  is  performed  to  form  the  template.  Then  the  template  is  modulated  near  a;„ 
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Spectrum  Square  Waveguide 


Frequency  GHz 

Figure  3-4:  Spectrum  of  field  in  dielectric  waveguide  using  two  dimensional  FDTD 
code  with  ky  =  600 

to  excite  the  simulation. 

To  validate  this  procedure  the  same  structure  of  Section  3.3  is  used.  The  FDTD 
grid  spacing  is  set  to  0.01  cm  and  the  time  step  set  to  the  Courant  limit  for  three 
dimensional  space  (2.45).  The  computational  domain  measured  37  x  41,  bounded 
with  a  first  order  Mur  absorbing  boundary  condition.  The  spatial  excitation  was  a 
two  dimensional  Gaussian  centered  inside  the  waveguide.  The  temporal  excitation 
was  a  narrow  Gaussian  pulse  with  =  16  as  described  in  Section  2.8.  The  excitation 
entered  the  simulation  volume  through  current  sources  at  the  desired  polarization  in 
this  case  Ez-  8192  time  steps  were  performed  and  the  last  4096  were  used  for  the 
Fourier  transform. 

The  spectrum  of  the  response  to  the  excitation  is  shown  in  Figure  3-4.  The  first 
peak  represents  the  temporal  frequency  of  the  fundamental  mode  that  corresponds 
to  the  pre-determined  propagation  constant  .  The  Fourier  transform  is  taken  at  the 
mode  frequency  to  produce  the  mode  template  in  Figure  3-5. 
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Figure  3-5;  Mode  template  for  Ef^  of  square  dielectric  waveguide  using  two  dimen¬ 
sional  FDTD  code. 

The  fundamental  mode  frequencies  at  different  propagation  constants  were  nu¬ 
merically  derived  with  the  FDTD  2D  code  and  plotted  against  the  analytic  curve  in 
Figure  3-6.  The  results  show  good  agreement  at  the  frequency  range  of  interest.  The 
deviation  at  the  higher  frequencies  is  due  to  the  numerical  dispersion  error  introduced 
by  the  discretization  of  Maxwell’s  equations  and  can  be  removed  with  a  smaller  grid 
size. 

3.5.2  Verification  of  Computational  Domain  Reduction 

Even  though  the  numerical  templates  lead  to  accurate  fundamental  mode  dispersion 
curves,  there  is  no  guarantee  that  the  computational  domain  can  be  reduced.  In  order 
to  show  this  reduction  several  FDTD  simulations  were  conducted  on  a  waveguide  with 
permittivity  e  =  2.8eo)  on  top  of  a  substrate  with  permittivity  e  =  2.2eo,  surrounded 
above  with  a  cladding  of  free  space.  The  substrate  is  mounted  on  a  ground  plane. 
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Analytic (d)  versus  Numerical { t) 


Figure  3-6:  Validation  of  mode  dispersion  curve.  Diamonds  represent  the  Marcatili 
calculation  and  the  triangles  represent  the  numerically  derived  results. 

The  basic  FDTD  simulation  volume  is  57  x  62  x  41  surrounded  by  an  8  layer  PML. 
The  grid  spacing  is  uniform  with  A  =  .001  m  and  the  time  step  is  set  at  the  3D 
Courant  limit.  A  modulated  Gaussian  pulse  with  /  =  20  GHz  and  /?  =  128  is  used 
with  either  a  numerical  template  or  a  2D  Gaussian  pulse  centered  in  the  waveguide 
with  widths  equal  to  the  width  and  height  of  the  waveguide.  The  center  point  is  set 
to  (18,15)  with  waveguide  dimensions  rh  =  lOA,  rw  =  lOA,  and  fh  =  0.  The  center 
of  the  perpendicular  section  of  the  waveguide  is  at  j  —  45  with  an  input  reference 
plane  at  j  =  22  and  an  output  reference  plane  at  z  =  41.  The  computational  domain 
is  varied  in  y  dimension  while  keeping  the  distances  between  the  reference  planes  and 
the  discontinuity  constant.  Two  larger  domains  of  57  x  82  x  41  and  57  x  92  x  41  were 
tested;  Figure  3-7  represents  a  top  view  of  the  configuration. 

Figure  3-8  and  Figure  3-9  represent  and  5'2i  calculated  at  each  of  the  three 
computational  domain  sizes  with  the  Gaussian  excitation.  A  detailed  description  of 
the  S  parameters,  and  5'2i,  can  be  found  in  Section  3.7.1.  Figure  3-10  and  Figure  3- 
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11  represent  and  S21  calculated  at  the  two  smallest  computational  domain  sizes 
using  a  numerical  mode  template  to  excite  the  problem.  The  extra  section  of  twenty 
grid  spaces  was  definitely  not  needed  when  using  the  mode  template  by  the  fact  that 
the  results  are  nearly  identical.  On  the  other  hand,  the  Gaussian  results  are  very 
different  between  the  small  and  middle  sized  domains.  Since  the  middle  and  large 
domains  are  nearly  identical,  the  mode  did  settle  to  an  acceptable  level  with  the 
addition  of  the  twenty  grid  spaces  in  the  y  dimension.  Thus,  using  a  mode  template 
in  this  simulation  reduced  the  computational  domain  by  over  20%.  In  what  would 
normally  take  two  hours,  the  computational  domain  reduction  means  a  CPU  savings 
of  over  20  minutes.  If  only  one  simulation  is  necessary,  since  finding  the  template 
takes  approximately  five  minutes,  the  savings  is  not  really  that  significant;  however, 
when  many  simulations  are  necessary  the  savings  can  be  substantial.  (All  simulations 
were  run  on  a  DEC3000  Server  Model  800  with  a  200  MHz  21064  Alpha  CPU.) 


3.6  Absorbing  Boundary  Condition  Demonstration 

Before  the  PML’s  introduction  and  its  subsequent  application  to  multi-layered  di¬ 
electrics,  a  first  order  Mur  ABC  was  usually  used  to  terminate  the  FDTD  computa¬ 
tional  domain  when  multiple  dielectric  media  axe  adjacent  to  the  domain  edges.  The 
second  order  Mur  cannot  be  used  because  of  discontinuities  at  dielectric/dielectric 
interfaces.  Figure  3-12  shows  the  time  domain  response  at  the  503rd  time  step  of  a 
simulation  to  analyze  the  effects  of  the  45°  bend  in  a  dielectric  rib  waveguide  with 
rw  =  10  mm,  rh  =  6  mm,  and  fh  =  4  mm.  The  cladding  is  free  space  and  the  per¬ 
mittivities  of  the  substrate  and  film  are  2.2eo  and  2.8eo  respectively.  In  open  space 
the  excitation  should  have  exited  the  simulation  volume  before  the  300th  time  step. 
Figure  3-12  clearly  shows  the  benefit  of  using  the  PML  over  the  first  order  Mur  ABC. 
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Figure  3-7:  Configuration  for  numerical  template  analysis 
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Figure  3-8:  Sn  as  measured  at  the  Input  Reference  Plane  for  an  analytical  Gaussian 
excitation  at  Source  Planes  A  (diamonds),  B  (triangles),  and  C  (stars). 


Figure  3-9:  ^21  as  measured  at  the  Output  Reference  Plane  for  an  analytical  Gaussian 
excitation  at  Source  Planes  A  (diamonds),  B  (triangles),  and  C  (stars). 
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Figure  3-10;  as  measured  at  the  Input  Reference  Plane  for  an  numeric  template 
excitation  at  Source  Planes  A  (diamonds),  and  B  (triangles). 


Figure  3-11;  S21  as  measured  at  the  Output  Reference  Plane  for  an  numeric  template 
excitation  at  Source  Planes  A  (diamonds),  and  B  (triangles). 
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Figure  3-12:  Time  domain  snap  shots  of  FDTD  simulation  with  1st  order  Mur  (left) 
and  PML  (right). 

3.7  Measurements 

Since  the  FDTD  is  a  time  domain  simulation  tool,  the  easiest  measurements  involve 
the  placement  of  software  probes  that  store  the  desired  field  quantities  at  each  time 
step.  When  frequency  measurements  are  needed,  a  sufficient  number  of  time  steps 
are  necessary  to  achieve  the  desired  frequency  resolution.  The  frequency  domain 
information  is  obtained  with  a  Fast  Fourier  Transform  (FFT)  of  the  time  domain 
data. 

3.7.1  Power  measurements 

To  analyze  waveguide  discontinuities,  the  S  parameters  are  found.  From  the  FDTD 
measurements  of  the  time-domain  waveforms,  the  S  parameters  may  be  calculated. 
The  first  step  is  to  insert  two  reference  planes  in  the  simulation  volume:  one  before 
the  discontinuity  and  the  other  after  the  discontinuity.  In  those  planes,  the  field 
components  are  saved  at  each  time  step.  Next,  the  Fourier  transform  of  each  field 
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component  at  each  location  is  taken  to  find  the  complex  time  harmonic  field  quan¬ 
tities.  Then  the  power  is  calculated  at  by  integrating  the  normal  component  of  the 
Poynting  vector  over  the  surface  of  the  reference  plane  and  taking  one  half  of  the  real 
part  [137]  such  that 


where  S  =  E  x  H*  ov 


and 


Power  =  J  ^  ^  H*ds'^  , 

(3.19) 

5,  =  EyHi  -  e,h;, 

(3.20) 

C  _  IP  TJ*  TP  IT* 

Uy  — 

(3.21) 

Q  _  TP  TJ*  771  TJ* 

(3.22) 

For  example,  in  the  dielectric  rib  waveguide  simulation,  the  propagation  is  in  the  y 
direction  so  the  power  transmitted  through  the  input  reference  plane  is  given  by 


Powevinput  plane  =  ^ AxA^Re  ^ E^H*  -  .  (3.23) 


3.7.2  S  Parameters 

The  scattering  parameters  can  be  derived  from  the  power  calculations.  However, 
first  the  incident  and  reflected  waves  must  be  separated.  This  is  done  by  running  two 
simulations;  one  with  a  straight  section  of  waveguide  with  no  discontinuity  (producing 
a  pure  incident  wave)  and  the  other  with  the  discontinuity  (producing  an  incident 
plus  a  reflected  wave).  The  waves  can  be  separated  by  subtracting  the  incident  wave 
calculated  in  the  first  simulation  from  the  incident  plus  reflected  waves  calculated  in 
the  second  simulation.  The  S  parameters  can  be  determined  from 


5n  =  101og(^) 


(3.24) 
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Figure  3-13:  Configuration  of  bends  tested,  from  left  to  right,  90°  bend,  90°  bend 
with  bevel,  and  45°  bend. 

and 

521  =  10  log  (3.25) 

where  Pj,  Pr,  and  Pt  represent  the  incident,  reflected  and  transmitted  power. 

3.8  Discontinuity  Analysis 

In  this  section  the  effects  of  various  types  of  bends  in  a  dielectric  rib  waveguide  are 
investigated.  The  different  bends  are  shown  in  Figure  3-13.  The  FDTD  simulation 
uses  a  uniform  grid  size  set  to  1  mm,  with  the  time  step  set  to  the  3D  Courant  limit. 
The  simulation’s  duration  is  1024  time  steps  and  the  subsequent  FFTs  use  the  data 
with  lx  padding. 

The  waveguide  is  constructed  with  a  film  layer  and  rib  of  permittivity  e  =  2.8eo) 
on  top  of  a  substrate  with  permittivity  e  =  2.2eo,  surrounded  above  with  a  cladding  of 
free  space.  The  substrate  is  mounted  on  a  ground  plane.  The  basic  FDTD  simulation 
volume  is  57  x  62  x  41  surrounded  by  an  8  layer  PML.  modulated  Gaussian  pulse 
with  /  =  18.8  GHz  and  (3  =  128  is  used  with  a  numerical  template  (Figure  3-14). 
The  center  point  is  set  to  (18,20)  with  waveguide  dimensions  rh  —  6A,  rw  =  lOA, 
and  fh  =  4 A.  The  center  of  the  perpendicular  section  of  waveguide  is  at  j  =  45  with 
an  input  reference  plane  at  j  =  22  and  an  output  reference  plane  at  z  =  41.  The 
loss  through  these  three  types  of  bends  can  be  seen  by  examining  ^21  in  Figure  3-15. 
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Figure  3-14:  Mode  template  for  Efi  < 
FDTD  code. 
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Figure  3-15:  ^21  rib  waveguide  through  three  different  types  of  turns:  90°  bend 
(diamond),  90°  bend  with  bevel  (triangle),  45°  bend  (star). 

None  of  the  turns  do  a  very  good  job  guiding  the  energy  around  the  turn.  The  45° 
bend  offers  the  best  performance  with  10  dB  of  loss.  One  of  the  advantages  of  using 
the  FDTD  method  in  this  type  of  analysis  is  that  snap  shots  can  be  taken  and  the 
progress  of  the  wave  can  be  visually  tracked.  Figures  3-16  and  3-17  show  how  the 
wave  travels  through  the  guide.  We  can  see  that  the  90°  bend  loses  most  of  its  energy 
in  radiated  out  the  back  of  the  turn  while  the  45°  bend  does  guide  energy  around  the 
turn  although  a  lot  is  still  radiated  out  through  the  back. 


3.9  Analysis  of  Rib  Geometry  on  Turn  Loss 

The  effects  of  the  film  height  to  rib  height  ratio  {fh  and  rf  in  Figure  3-1)  are  inves¬ 
tigated  in  this  section.  Figures  3-18  through  3-22  show  the  mode  templates  for  the 
five  guides  tested.  Each  guide  is  excited  with  modulated  Gaussian  pulse.  The  guides 
were  10  mm  wide  and  a  45°  bend  was  used  as  the  discontinuity  tested.  The  results 
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Figure  3-16;  Time  domain  snap  shots  of  propagation  through  a  90°  bend  at  173,  233, 
293,  and  323  time  steps.  (Clockwise  starting  at  upper  left.) 
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Figure  3-17:  Time  domain  snap  shots  of  propagation  through  a  45°  bend  at  173,  233, 
293,  and  323  time  steps.  (Clockwise  starting  at  upper  left.) 
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Figure  3-18:  Mode  template  for  E^  ^  of  dielectric  rib  waveguide  with  film  height: 
6  mm,  rib  height:  4  mm  and  rib  width:  10  mm. 


Figure  3-19:  Mode  template  for  E^i  of  dielectric  rib  waveguide  with  film  height: 
5  mm,  rib  height:  5  mm  and  rib  width:  10  mm. 


show  that  the  smaller  film  height  to  rib  height  ratio  the  better  the  bend  will  guide 
the  energy.  Over  10  dB  of  improvement  can  be  seen  when  the  ratio  is  decreased  from 
1.5  to  0.  In  other  words,  a  rectangular  waveguide  on  top  of  a  substrate  has  better 
performance  than  a  rib  guide.  And,  as  the  film  thickness  is  reduced  compared  to 
the  film  height  the  better  the  performance.  Better  performance  is  attributed  to  the 
open  nature  of  the  film  layer.  As  the  rib  height  decreases  relative  to  the  film  height, 
the  guide  approaches  a  dielectric  slab  waveguide  that  would  offer  no  guidance  in  the 
transverse  direction. 
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Figure  3-20:  Mode  template  for  Ef  ^  of  dielectric  rib  waveguide  with  film  height: 
4  mm,  rib  height:  6  mm  and  rib  width:  10  mm. 


Figure  3-21:  Mode  template  for  El^  of  dielectric  rib  waveguide  with  film  height: 
3  mm,  rib  height:  7  mm  and  rib  width:  10  mm. 


Figure  3-22:  Mode  template  for  Ef  ^  of  dielectric  rib  waveguide  with  film  height: 
0  mm,  rib  height:  10  mm  and  rib  width:  10  mm. 
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Figure  3-23:  521  through  45°  bend  of  dielectric  rib  waveguide.  Starting  from  top 
curve  to  bottom  curve  the  film  height  to  rib  height  ratios  are  0.0,  0.43,  0.67,  1.0,  and 
1.5. 

3.10  Discussion 

The  Finite  Difference  Time  Domain  numerical  technique  is  an  excellent  tool  to  study 
the  dielectric  rib  waveguide.  The  openness  of  the  structure  can  be  simulated  with 
the  Berenger  Perfectly  Matched  Layer  if  each  dielectric  layer  is  matched  considering 
all  the  layers  within  the  simulation  volume.  The  computation  domain  can  be  reduce 
if  the  fundamental  mode’s  spatial  distribution  is  calculated  first  and  used  to  excite 
the  mode.  The  better  excitation  means  that  shorter  distances  between  the  excitation 
plane  and  the  discontinuity  are  needed  for  the  mode  to  settle.  A  twenty  percent 
reduction  was  shown  in  Section  3.5.2.  This  reduction  implies  bigger  problems  can  be 
simulated  within  the  same  time  as  those  simulations  using  other  excitation  techniques, 
or  the  same  simulations  can  be  run  in  shorter  times. 

Although  the  dielectric  rib  waveguide  provides  very  low  loss  for  high  frequency 
signals  compared  to  microstrip  or  coplanar  structures,  the  rib  waveguide  bend  does 
introduce  a  significant  amount  of  loss.  As  the  abruptness  on  the  turn  is  lessened,  the 
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rib  guides  more  energy  around  turn;  however,  gradual  turns  require  more  space  on 
an  integrated  circuit.  With  the  improvements  to  the  FDTD  method  presented  in  this 
chapter,  rib  bend  structures  can  be  studied  to  improve  their  loss  characteristics. 
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Chapter  4 


Impedance  Boundary  Condition 
for  Thin  Finite  Conducting  Sheets 
in  FDTD 

4.1  Introduction 

In  FDTD  simulations  involving  highly  conductive  materials,  such  as  the  metal  case  of 
a  computer  system,  the  tangential  electric  fields  are  typically  set  to  zero  on  the  surface 
of  the  material.  This  Perfect  Electric  Conductor  (PEC)  assumption  ignores  any  loss 
associated  with  the  less  than  infinite  conductivity.  The  errors  that  result  from  this 
approximation  are  considered  against  the  large  cost  of  discretizing  the  lossy  material. 
Since  the  wavelength  of  a  highly  conductive  material  is  very  small,  in  order  to  capture 
the  loss  mechanism  within  the  conductor  the  simulation  volume  must  be  divided  into 
a  very  fine  grid.  When  the  brute  force  method  is  used,  computer  resources  are  wasted 
by  having  to  finely  divide  the  the  other  less  dense  materials  {e.g.  free  space).  In  these 
materials  the  fine  grid  is  not  necessary  to  capture  the  physics  of  the  problem.  The 
computational  size  can  be  reduced  with  a  sub-gridding  scheme  where  the  conducting 
material  uses  much  smaller  grid  than  the  rest  of  the  computational  domain  [9,  10,  12]. 


107 


108  CHAPTER  4.  IBC  FOR  THIN  FINITE  CONDUCTING  SHEETS  IN  FDTD 


For  very  good  conductors;  however,  the  grid  size  is  so  small  that  even  sub-gridding 
is  not  a  viable  option.  To  overcome  the  resource  problem,  the  surface  of  the  highly 
conductive  material  can  be  replaced  with  an  Impedance  Boundary  Condition  (IBC) 
[126]-[130].  An  IBC  is  only  appropriate  when  the  simulation  volume  of  interest  is  on 
one  side  of  the  conductive  material.  However,  IBCs  have  the  added  complication  that 
they  are  usually  frequency  dependent  and  are  not  directly  applicable  to  the  the  stan¬ 
dard  frequency  independent  FDTD  equations.  In  this  case  the  FDTD  equations  must 
be  modified  to  incorporate  the  dispersive  nature  of  the  surface  [109]-[120].  Typical 
frequency  domain  equations  relating  the  electric  fields  and  magnetic  fields  become 
convolution  equations  is  time.  Convolutions  have  a  large  computational  and  memory 
overhead,  but  if  the  convolution  integral  can  be  approximated  by  a  sum  of  exponen¬ 
tials  then  recursive  convolution  can  be  implemented  and  the  memory  requirements 
reduced. 

A  synthetic  conductivity  has  been  used  with  normal  FDTD  difference  equations 
for  lossy  media  to  model  good  conductors  [117].  The  synthetic  conductivity  is  derived 
by  comparing  the  numerical  impedance  of  the  difference  equations  and  the  actual  im¬ 
pedance  of  a  good  conductor  at  specified  frequency.  In  this  way  the  derived  synthetic 
conductivity  is  inserted  into  the  FDTD  equations  at  the  boundary  only;  thus  being  a 
surface  boundary  condition  similar  to  an  IBC.  The  advantage  to  this  method  is  that 
no  new  equations  are  needed;  however,  it’s  major  disadvantage  is  that  it  is  good  at  a 
single  frequency. 

An  Impedance  Boundary  Condition  that  incorporates  the  frequency  dependence 
was  developed  for  thin  dielectric  coatings  over  PEC  surfaces  [127].  The  IBC  starts 
with  the  analytically  derived  expression  for  the  impedance  in  the  frequency  domain. 
The  expression  is  expanded  in  a  Taylor  series  about  the  wave  number,  k,  and  trans¬ 
formed  to  the  time  domain.  The  resulting  partial  differential  equation  is  fourth  order 
in  time  and  first  order  in  space.  The  difference  equations  do  model  the  frequency 
characteristics  fairly  well;  however,  the  equations  are  quite  complicated  and  rely  on 
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many  past  values  of  the  surface  electric  field.  Furthermore,  the  condition  is  good  at 
normal  incidence  only  thus  limited  to  ID  structures.  This  work  is  extended  to  in¬ 
clude  non-normal  incidence  [129] ;  however,  the  modified  IBC  is  only  good  for  a  single 
incidence  angle. 

Like  the  IBC  of  [127],  Tassoudji  [128]  starts  with  a  frequency  domain  expression 
for  the  impedance.  This  time  for  a  thin  sheet  of  a  very  good  conductor.  The  frequency 
domain  equation  is  transformed  to  the  Laplace  domain  and  the  subsequent  expression 
is  approximated  by  a  sum  of  first-order  rational  functions.  This  approximation  is  crit¬ 
ical  to  the  IBC’s  usefulness  in  FDTD  simulations.  The  first-order  rational  functions 
are  exponentials  in  time  and  thus  a  recursive  convolution  scheme  is  implemented. 
This  FDTD  IBC  scheme  has  been  applied  to  one-dimensional  (ID)  FDTD  simula¬ 
tion  of  a  plane  wave  reflecting  from  35 //m  thick  copper  plate  with  conductivity  of 
cr  =  5.8  X  10^  Compared  with  the  anal3fi;ic  solution,  the  results  have  shown 

excellent  agreement  up  to  30  GHz.  The  benefit  of  this  scheme  is  that  it  is  good  over 
a  wide  range  of  frequencies. 

The  same  approach  was  applied  to  lossy  dielectric  portions  of  a  FDTD  computa¬ 
tional  domain  with  an  IBC  where  the  frequency  domain  impedance  expression  was 
approximated  by  a  rational  function  using  Mathematica  and  then  expanded  using 
partial  fraction  expansion.  In  this  way,  like  [128]  a  recursive  convolution  can  be  used. 

In  this  chapter  the  IBC  for  thin  finite  conducting  sheets  of  [128]  is  examined  and 
extended  to  two  dimensions.  In  order  for  this  IBC  to  be  a  practical  FDTD  tool,  this 
extension  is  very  important.  Also,  guidelines  for  the  number  of  expansion  terms  are 
needed  in  order  to  give  users  an  a  priori  estimate  of  the  resources  required  to  perform 
a  given  simulation. 

This  chapter  describes  the  ID  Impedance  Boundary  Condition  for  thin  finite  con¬ 
ducting  sheets  and  recursive  convolution.  Then  the  accuracy  of  the  ID  IBC  is  exam¬ 
ined  and  quantified.  Finally,  the  IBC  is  extended  to  two  dimensions  and  validated 
through  the  measurements  of  the  quality  factor  of  rectangular  resonant  cavities  at 
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Figure  4-1:  ID  transmission  line  model  used  for  IBC  formulation. 


different  frequencies  and  conductivity  values. 


4.2  Derivation  of  Impedance  Boundary  Condition 

To  construct  an  impedance  boundary  condition  we  start  with  the  impedance  equation 
that  relates  the  tangential  components  of  the  electric  and  magnetic  field  at  an  interface 
between  two  different  media, 


Et&xiiN)  ~  E{oj) 


(4.1) 


where  Stan  and  H^an  are  the  tangential  components  of  the  electric  and  magnetic  fields 
and  n  is  the  normal  unit  vector  pointing  out  at  the  interface  [128].  For  a  conducting 
media  with  properties  e,  fXg  and  a  representing  the  permittivity,  permeability,  and 
conductivity  respectively,  if  the  media  is  a  good  conductor  such  that  1  «  -^  then 
the  permittivity  of  the  conductor  can  be  approximated  by 


e 


a 

i— 

U) 


(4.2) 


In  one  dimensional  simulations,  since  propagation  is  normal  to  the  surface,  the 
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surface  impedance  of  a  thin  good  conductor  can  be  modeled  as  a  transmission  line 
with  thickness  I  and  an  intrinsic  impedance  rjc  =  Prom  transmission  line  theory, 
if  a  transmission  line  is  terminated  with  impedance  rit  then  the  impedance  seen  at 


distance  I  from  the  termination  is 


^in  Vc 


rjt  -  r)cteinh{ikcl) 
rjc  -  rjt  tanh(zA:cO  ’ 


When  the  termination  impedance  is  free  space  i.e.  r]t  =  t]o  ^  r]c  ,  then  input  impe¬ 


dance  can  be  approximated  by 


Zin  ~  r)cCoth{ikcl). 


In  order  to  utilize  this  impedance  condition,  it  is  necessary  to  transform  equation  into 
the  Laplace  domain.  The  Laplace  transform  of  (4.4)  is 


Zinis)  =  y^coth  ■ 


Using  power  series  expansions  for  cosh  and  sinh  [141], 


cosh  ^  =  JJ  1  + 


{2k  —  1)%" 


sinh  z  =  zY[  1  + 


the  input  impedance  (4.5)  can  be  represented  by  quotient  of  two  infinite  products: 


Zin{s^  j 


(2m— l)27r2 


n  [i  _| _ inooJL- 

1  Jii  U  ^  (2m-l)^7r^ 

n  [i  +  ^s] 

m—1 
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To  use  this  relationship  within  the  FDTD  framework  the  impedance  condition  (4.3) 
must  be  transformed  into  the  time  domain, 


Et^n{t)  =  J  E{t-  t)  [n  X  JTtan(T) 


dr. 


(4.9) 


A  simple  expression  for  .Etan(i)  niust  be  found  that  can  be  converted  to  update  equa¬ 
tions  that  do  not  explicitly  use  the  convolution  integral  of  (4.9).  This  can  be  accom¬ 
plished  by  approximating  (4.8)  with  P  terms  in  the  denominator  and  P  -1  terms 
in  the  numerator.  The  quotient  can  be  expressed  as  a  sum  of  quotients  using  partial 
fraction  expansion  and  then  transformed  into  the  time  domain.  The  result  is  a  sum 
of  exponentials, 

Z{t)  ^  f;  Cme^-\  (4.10) 

m—l 

where  is  the  mth  pole  and  Cm  in  the  mth  residue  of  (4.8),  given  by 


Am 


li'oCFP 


(4.11) 


and 


a 


Zin{s') 


S — A.m 


(4.12) 


4.3  Recursive  Convolution 


Recursive  convolution  is  possible  when  one  of  the  operands  is  an  exponential  function. 
Since  we  have  expressed  the  impedance  of  (4.9)  at  a  sum  of  exponentials  recursive 
convolution  can  be  used.  As  in  [130],  we  begin  with  the  following  convolution  integral. 
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and  discretize  it  so  that 

n-l  f{k+l)At  ,,  ^ 

y{nAt)  =  X!  /  ^^a;(r)dr. 

A;=o 

Assuming  that  x{t)  is  constant  over  the  interval  {i.e.  piecewise  constant), 
gration  over  r  yields 

y{nAt)  =  -  1). 

^  fc=0 

Remove  the  nth  term  from  the  summation, 
y{nAt)  = 

^  k=0 

+  ^x{kAt)e-^^\e’’^^  -  1), 

0 


and  rewrite  the  sum  over  n  —  2  so  that 

-  1) 

fc=0 

=  -  1) 

k=0 

=  -  1) 

k=0 

=  e-“‘!/((n  -  1)A(). 


The  result  is  the  basis  for  a  recursive  convolution  equation: 

y(nAt)  =  e"^^*?/((n  -  l)At) 

+  ^x((n-l)At)(l-e-'’^*). 
0 


(4.14) 
the  inte- 

(4.15) 


(4.16) 


(4.17) 


(4,18) 


In  the  derivation  above  we  assume  x{t)  could  be  approximated  as  x{nAt)  over  the 
interval  [nAt,  (n  +  l)At].  If  we  assume  a  piecewise  linear  approximation  instead  of 
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piecewise  constant,  x{t)  is  approximated  as 

x{t)  ~  x{nAt)  +  +  1)^^)  “  x{nAt))  (4.19) 


so  that 


f{k+l)At 

y{nAt)  =  XI  / 

k=Q 


and  the  recursion  formula  becomes 


^  _  Ic 

x{nAt)  H - — — (x((n  +  l)At)  —  x{nAt)) 


At 


dr 
(4.20) 


y{nAt)  =  e  ^^*y{{n  —  l)At) 

+  f  {x(nAt)  [l  +  ^(e-**^*  -  1)]  (4.21) 

+  x((n  -  1) At)  -  e-'’^‘(l  +  ^)] }. 


4.4  Implementation  of  Impedance  Boundary  Con¬ 
dition  in  FDTD 


Now  that  we  have  the  form  of  the  recursive  convolution  integrals  we  can  construct 
FDTD  IBC  update  equations. 


4.4.1  Impedance  Equations  -  Piecewise  Constant  Assump¬ 
tion 

The  derivation  of  the  previous  section  can  be  used  to  construct  the  FDTD  update 
equations  if  we  let 


y[nAt)  ^  GZ, 
a;((n-l)At)  n  x  S^nK 


(4.22) 

(4.23) 
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and 


^  AmnAt 


(4.24) 


Substitution  into  (4.19)  gives 

G«  =  (4.25) 

+  (ft  X  -  1)  (4-26) 

The  tangential  electric  field  update  equation  at  the  conductor  surface  is 

A".»  “Eg;  (4.27) 

m=l 

The  impedance  equation  relates  the  electric  and  magnetic  fields  at  the  same  position 
and  at  the  same  time.  However,  in  this  work,  we  use  the  tangential  electric  field  on 
the  surface  accompanied  by  the  tangential  magnetic  field  half  a  grid  space  away  and 
we  take  a  time  average  of  the  electric  field  so  that  the  update  equation  becomes 

^»=-.^;‘+2l;G»t  (4.28) 

m=l 

The  tangential  electric  field  on  the  surface  at  the  current  time  step,  n,  is  a  function  of 
last  tangential  electric  field  and  the  last  sum,  Gm-  Compared  to  the  full  convolution, 
memory  requirements  are  low  since  this  method  requires  only  2P  variables  for  the 
poles  and  residues  and  P  variables  for  each  grid  point  on  the  surface  containing  a 
tangential  electric  field. 
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4.4.2  Impedance  Equations  -  Piecewise  Linear  Assumption 

When  more  accuracy  is  desired  the  piecewise  linear  assumption  can  be  used  so  that 
the  tangential  electric  fields  update  equation  is  still  (4.28)  but  now 


G 


n 

m 


+ 

+ 


^n-l  Am^t 

(ft  X  H^J)^ 

■  1 

Am  At 
-  1 

(ft  X 

.AmAt 

-(e^™A*  _  1)  _  1 


+  e' 


‘(1 


Am  At' 


(4.29) 


Improved  accuracy  is  obtained  at  the  expense  of  a  more  complicated  equation  and 
an  additional  memory  location  for  each  tangential  electric  field  component  on  the 
surface. 


4.5  Extension  of  IBC  to  Two  Dimensions 

The  use  of  the  transmission  line  to  model  the  impedance  of  a  thin  finite  conducting 
sheet  is  completely  justified  in  ID  simulations  since  propagation  is  fixed  in  the  normal 
direction.  On  the  other  hand,  since  the  the  angle  of  incidence  is  not  restricted  to 
normal  incidence  in  2D  simulations,  we  must  verify  that  the  model  is  still  valid. 

Consider  phase  matching  at  the  surface  of  the  conductor,  located  at  2:  =  0,  where 
ki  is  the  wave  number  on  the  free  space  side  with  incidence  angle  6i  and  kc  is  the  wave 
number  in  the  conductor  with  a  transmission  angle  of  dc-  The  incident  components 
of  the  k  vector  are  kix  =  ko  sin  6i  and  ki^  =  ko  cos  9i  and  the  conductor  wave  number 
is  kc  =  ko^l  +  i±-^. 

Since  phase  matching  says  kix  =  kcx,  we  know  k^^  =  kl-  kf^.  Solving  for  kcz, 


kcz  =  ^0  fcos6>i^  +  )  , 


(4.30) 
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and  using  the  good  conductor  approximation,  we  find 

The  angle  in  the  conductor  is  given  by 

tan^c  =  p  jT  -.  =  sing; (4.32) 

Re{A;c4  V 


Since  is  very  large,  the  propagation  direction  will,  for  all  practical  purposes, 
be  normal  to  the  surface  regardless  of  the  incidence  angle.  Therefore,  we  expect  the 
transmission  line  model  to  apply  even  in  the  two  dimensional  case. 

Figure  4-2  is  a  schematic  of  a  2D  rectangular  resonator  with  finite  conducting 
walls.  The  relationships  between  the  tangential  electric  and  magnetic  walls  and  the 
surface  normal  are  shown  and  are  used  to  construct  IBC  equations  for  each  of  the 
four  walls.  If  the  computational  domain  is  Nx^  N^,  the  piecewise  constant  equations 
are 


=  -E^-\0,k) 

+  2  f;  -  1) 

m=l  ^  ^ 


(4.33) 


for  the  left  wall; 


E^{Nx,k)  =  -E^-\Nx,k) 


+ 


y 
p 

2E 

m— 1 


(4.34) 
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n-l  AmAt 
m  ^ 


Cr, 


+  HrHNx-ik)^{e 


Am^t 


1) 


for  the  right  wall; 


(4.35) 
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Figure  4-2:  2D  configuration  used  for  IBC  equations  of  a  rectangular  resonator. 


+  2^ 


-1)  - 


for  the  lower  wall;  and 


E^{i,N,)  =  -E^-\i,N,)  (4.36) 

+  2  X;  -  HTHi,N,  -  i,  -  1) 

m=l  ^  ^ra 

for  the  upper  wall. 

The  piece  linear  equations  can  easily  be  constructed  by  modifying  equations  (4.34) 
through  (4.37)  keeping  equation  (4.30)  in  mind. 


4.6  Effects  of  FDTD  Simulation  Parameters  on  IBC 
Accuracy 


In  this  section  a  numerical  study  of  the  effects  of  the  FDTD  simulation  parameters 
such  as  grid  spacing,  time  step,  and  the  number  of  expansions  terms  is  presented. 
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Figure  4-3:  ID  configuration  for  IBC  error  analysis. 


4.6.1  ID  FDTD  Simulation  Configuration 


Figure  4-3  is  a  representation  of  the  ID  FDTD  simulation  set  up  used  to  analyze  the 
effects  of  FDTD  parameters  on  the  IBC  accmacy.  An  unmodulated  Gaussian  pulse  is 
used  to  excite  the  problem  at  one  end  of  the  computational  domain,  point  A,  with  the 
boundary  condition  E^{k  =  0)  =  g{n).  The  IBC  is  placed  at  point  B  and  a  first  order 
Mur  boundary  condition  is  placed  at  the  other  end  of  the  computational  domain,  point 
C.  The  propagation  direction  is  2  and  a  TEM  wave,  {Ex,  Hy)  is  launched  toward  the 
IBC.  The  computational  domain  is  3000A  long  with  B  in  the  middle  at  z  =  1500A 
or  more  simply  k  =  1500.  As  described  in  Chapter  2,  the  time  step  is  set  to 


where  c  is  the  speed  of  light  and  CFL  a  real  number  greater  than  or  equal  to  the 
Courant  limit.  The  normal  ID  FDTD  algorithm  requires  CFL  >  1  for  stability.  Two 
simulations  are  run  for  any  given  test  set  up.  The  first  simulation  does  not  include 
the  IBC  at  A:  =  1500.  Instead  the  field  is  allowed  to  freely  propagate  to  the  right 
and  the  field  is  measured  at  A:  =  1499  and  saved  at  all  time  steps  to  be  used  as 
the  incident  field  when  calculating  the  reflection  coefficient  of  the  IBC.  The  second 
simulation  also  measures  the  electric  field  at  k  =  1499;  however  this  time,  the  IBC  is 
in  place  at  A:  =  1500.  Here  the  stored  field  includes  both  the  incident  and  reflected 
waves.  By  using  this  total  field  and  the  stored  incident  field  the  reflection  coefficient  is 
calculated  up  to  10  GHz.  Even  though  the  reflection  coefficient  would  be  a  sufficient 
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Figure  4-4:  Transmissivity  FDTD  (dashed)  vs.  Analytic  (solid)  for  a  35 /xm  thick 
sheet  of  copper  with  cr  =  5.8  x  P  =  20,  A  =  .005  m,  CFL  =  2.0,  and 

^  =  50. 


measure  of  the  IBC,  the  transmissivity  is  calculated  from  the  reflection  coefficient  via 

t(/)  =  101og[l-|P(/)n.  (4.37) 

The  transmissivity  has  the  nice  property  that  it  increases  with  increasing  frequency 
in  this  application. 

The  frequency  content  of  the  incident  pulse  is  controlled  through  the  selection 
of  /?  from  Chapter  2  once  the  time  step  is  determined  by  A  and  CFL.  Figure  4-4 
is  a  plot  of  the  transmissivity  measured  from  a  35  iim  thick  sheet  of  copper  with 
(7  =  5.8  X  10^r2“^m"'h  The  IBC  is  within  IdB  up  to  3  GHz  with  only  20  expansion 
terms.  In  the  next  section  the  FDTD  parameters  are  varied  and  the  effects  studied. 
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4.6.2  Transmissivity  Error  Measurements 


The  error  analysis  is  based  on  a  direct  comparison  between  the  FDTD  transmissivity 
results  versus  the  analytical  transmissivity  calculation  using 


e  = 


Eill[4(/i)  -  tFDTpjfiW 
N 


(4.38) 


where  tpoToifi)  is  the  transmissivity  calculated  from  the  FDTD  simulation  and  ta{fi) 
is  the  equivalent  analytical  calculation  at  frequency,  /j. 

The  first  FDTD  parameter  examined  is  the  time  step.  In  [128],  the  time  step  was 
fixed  relative  to  the  grid  spacing  with  CFL  =  2.0.  However,  the  FDTD  stability  cri¬ 
teria  places  the  lower  limit  at  CFL  =  1.0.  A  larger  time  step  can  reduce  the  run  time 
and  improve  the  frequency  resolution  in  any  post  simulation  processing.  Typically, 
FDTD  codes  allow  the  users  to  set  the  grid  spacing  while  the  code  automatically  sets 
the  time  step  at  the  Courant  limit,  so  it’s  important  to  know  if  the  IBC  requires 
a  time  step  smaller  than  the  Courant  stability  limit  in  order  to  bound  the  error. 
Figure  4-5  shows  the  results  from  a  numerical  experiment  with  two  very  different 
simulation  configurations.  The  conductivities  are  different,  5.8  x  10^  versus 

5.8  X  10^  the  number  of  expansion  terms  are  different,  P  =  3  versus  P  =  125; 

and  finally  the  grid  spacing  is  different,  A  =  .005  m  versus  A  =  .00125  m.  Both  sim¬ 
ulations  used  a  sheet  thickness  of  35  /xm  and  the  (3  is  adjusted  at  each  simulation  to 
ensure  the  same  frequency  content  in  the  excitation.  Also,  all  simulations  were  run 
for  8192  times  steps  and  the  Fourier  transforms  were  taken  with  padding  at  16  times 
the  original  length.  The  error  is  calculated  with  (4.38).  The  results  show  that  in 
order  to  keep  the  error  low  the  IBC  definitely  has  a  stricter  time  step  constraint  than 
the  FDTD  Courant  stability  limit.  In  fact,  both  ID  experimental  configurations  had 
the  same  CFL  limit.  The  results  show  that  CFL  >  1.82  is  necessary  to  bound  the 
error.  This  does  not  mean  the  algorithm  is  unstable  when  CFL  is  set  to  the  Courant 
limit  it  just  means  the  accuracy  is  very  poor.  Furthermore,  it  is  important  to  note 
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CFL 

Figure  4-5:  Error  vs.  CFL  for  a  35  fj,m  thick  conductor  for  two  different  FDTD 
configurations:  (1)  cr  =  5.8  x  P  =  3,  A  =  .005  m  (diamond)  and  (2) 

cr  =  5.8  X  10^O“^m~^,  P  =  125,  A  =  .00125m  (triangle). 

that  the  CFL  is  even  larger  than  the  3D  FDTD  Courant  limit  of  ^/3  for  a  uniform 
grid. 

Next  the  effects  of  grid  size  and  number  of  expansion  terms  are  studied.  Figures  4- 
6  through  4-11  are  plots  of  the  transmissivity  error  of  equation  (4.38)  versus  the 
number  of  expansion  terms  at  three  different  grid  spacings,  A  =  .005  m,  A/2,  and 
A/4.  All  simulations  used  a  conductor  thickness  of  35 //m  and  the  P’s  were  50,  100, 
and  200  respectively.  Figure  4-9  are  the  results  with  copper,  a  —  5.8  x  10^,  while 
Figure  4-6,  Figure  4-7,  and  Figure  4-8  are  the  results  for  a  conductivities  three  orders 
of  magnitude  less  than  copper  and  Figure  4-10  and  Figure  4-11  are  the  results  for  a 
conductivities  2  orders  of  magnitude  greater  than  copper. 

The  results  show  the  general  trend  that  smaller  grid  size  will  reduce  the  error. 
This,  of  course,  makes  sense.  The  FDTD  algorithm  is  more  accurate  with  a  smaller 
grid  size.  Also,  the  impedance  in  the  FDTD  IBC  is  approximated  by  electric  and 
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magnetic  fields  that  are  half  a  grid  space  apart,  so  the  closer  the  grid  space  the  more 
accurately  the  IBC  scheme  will  approximate  the  actual  impedance.  However,  the 
numerical  results  also  show  that  the  number  of  expansion  terms  plays  a  big  role  in 
the  accuracy  of  the  IBC  algorithm  and  more  terms  does  not  necessarily  mean  more 
accuracy.  For  example,  in  Figure  4-6  we  see  that  the  error  levels  out  at  about  20 
expansion  terms.  Interestingly,  the  optimum  setting  for  P  is  about  5.  Although  the 
differences  in  error  are  less  than  .2dB  and  .IdB  for  A  =  .005m  and  A  =  .00125m 
respectively,  these  errors  would  accumulate  in  simulations  where  multiple  reflections 
from  the  surface  occur. 

The  optimum  number  of  expansion  terms  is  a  function  of  conductivity  and  grid 
spacing.  As  Figure  4-9  indicates,  the  best  values  for  P  are  approximately  50,  80,  and 
120  for  a  conductivity  of  a  =  5.8  x  10^  and  a  grid  spacing  of  A  =  .005  m,  A  =  .0025  m, 
and  A  =  .00125  m  respectively.  Whereas  for  low  and  high  conductivities  the  effect  of 
grid  spacing  is  much  reduced. 

The  results  clearly  show  that  as  the  conductivity  increases  the  number  of  terms 
need  to  accurately  approximate  the  impedance  increases.  In  fact  for  the  largest 
conductivity  shown,  cr  —  5.8  x  10®,  the  IBC  requires  more  than  250  terms  to  keep  the 
error  below  1  dB. 

4.6.3  Optimum  Settings 

The  fact  that  there  is  an  optimum  number  of  terms  for  the  piecewise  constant  IBC 
seems  to  go  against  intuition.  If  additional  terms  more  accurately  model  the  impe¬ 
dance  then  why  doesn’t  the  error  decrease  as  the  number  of  expansion  terms  increase? 

We  begin  with  a  look  at  the  error  results  for  the  piecewise  linear  IBC  for  copper 
with  cr  =  5.8  X  10’’  n“’m“’  at  three  different  grid  spacings,  as  shown  in  Figure  4-12 
versus  the  piecewise  constant  error  in  Figure  4-13.  The  two  plots  are  combined  in 
Figure  4-14.  The  results  show  that  there  is  no  optimum  number  of  expansion  terms 
(at  least  up  to  300  terms)  when  using  the  piecewise  linear  IBC.  Next,  we  look  at  the 
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Figure  4-6:  Error  vs.  expansion  terms  for  a  35  fim  thick  conductor  with  cr  =  5.8  x 
A  =  .005m  (diamond),  A  =  .0025m  (triangle),  A  =  .00125m  (star). 


Figure  4-7:  Error  vs.  expansion  terms  for  a  35 /im  thick  conductor  with  cr  =  5.8  x 
10®J7~^m“^;  A  =  .005m  (diamond),  A  =  .0025m  (triangle),  A  =  .00125m  (star). 
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Figure  4-8:  Error  vs.  expansion  terms  for  a  Zbiim  thick  conductor  with  cr  =  5.8  x 
A  =  .005m  (diamond),  A  =  .0025m  (triangle),  A  =  .00125m  (star). 


Figure  4-9:  Error  vs.  expansion  terms  for  a  35  (xm  thick  conductor  with  cr  =  5.8  x 
A  =  .005m  (diamond),  A  =  .0025m  (triangle),  A  =  .00125m  (star). 
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Figure  4-10:  Error  vs.  expansion  terms  for  a  35  jim  thick  conductor  with  <7  =  5.8  x 
10^  A  =  .005m  (diamond),  A  =  .0025m  (triangle),  A  =  .00125m  (star). 


Figure  4-11:  Error  vs.  expansion  terms  for  a  35 /xm  thick  conductor  with  cr  =  5.8  x 
10^f^“^m“^;  A  =  .005m  (diamond),  A  =  .0025m  (triangle),  A  =  .00125m  (star). 
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magnitude  of  main  term  in  the  piecewise  constant  IBC  formulation,  —  1,  and 

plot  it  versus  expansion  number  at  the  three  time  steps  of  Figure  4-12  as  shown  in 
Figure  4-15.  Figure  4-15  shows  that  we  should  expect  a  leveling  off  of  the  error  since 
the  IBC  term  levels  off  rather  quickly.  Since  the  exponential  term  is  proportional  to 
the  square  of  the  expansion  term  number,  it  follows  that  the  influence  of  the  term 
will  quickly  decrease.  Furthermore,  the  other  factors  in  the  exponential  term  should 
show  similar  results.  As  the  conductivity  decreases  the  number  of  strong  influential 
terms  decreases  (Figures  4-6  through  4-11).  As  the  conductor  thickness  increases  the 
number  of  influential  terms  also  increases  (Figure  4-16). 

There  are  essentially  two  approximations  that  contribute  to  the  transmissivity  er¬ 
ror.  The  first  is  from  the  approximation  of  the  impedance  by  a  sum  of  exponentials, 
and  the  other  is  from  the  approximation  to  the  convolution  integral.  The  two  con¬ 
volution  approximations  used  are  piecewise  constant  and  piecewise  linear.  When  few 
expansion  terms  are  used,  the  impedance  error  is  the  dominant  error.  As  the  number 
of  terms  increases,  the  impedance  error  decreases.  If  fact,  it  should  approach  zero  as 
m  oo.  Therefore,  as  m  becomes  sufficiently  large,  the  convolution  integral  approx¬ 
imation  error  will  dominate  the  overall  error  and  increasing  m  will  not  improve  the 
error.  This  is  clearly  evident  in  Figure  4-14  as  seen  by  the  leveling  of  the  error  curves 
at  some  small  but  finite  error  as  m  increases.  The  optimum  settings  evident  in  the 
piecewise  constant  must  result  from  the  combination  of  the  two  types  of  errors,  and 
the  piecewise  linear  error  is  evidently  sufficiently  small  such  that  it  does  not  produce 
an  optimum  setting. 


Depending  on  the  problem  being  considered,  a  non-optimum  number  of  expansion 
terms  may  not  lead  to  significant  error  since  the  errors  are  still  small.  However,  if  the 
problem  includes  many  reflections  from  the  IBC  as  in  resonant  cavities,  the  optimum 
number  of  terms  should  be  considered  when  using  the  piecewise  constant  formulation. 
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Figure  4-12:  Piecewise  linear  error  vs.  expansion  terms  for  a  35  /rm  thick  conductor 
with  (T  =  5.8  X  llT  A  =  .005  m  (solid),  A  =  .0025  m  (dash),  A  =  .00125  m 

(dash-dot). 


Figure  4-13:  Piecewise  constant  error  vs.  expansion  terms  for  a  35  ixm  thick  conductor 
with  cr  =  5.8  X  10^fi"^m“^;  A  =  .005m  (solid),  A  =  .0025m  (dash),  A  =  .00125m 
(dash-dot). 
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Figure  4-14:  Comparison  of  error  behavior  of  piecewise  constant  IBC  and  piecewise 
linear  IBC  for  a  35  yum  thick  conductor  with  cr  =  5.8  x  10^  A  =  .005  m 

(solid),  A  =  .0025m  (dash),  A  =  .00125m  (dash-dot). 


Figure  4-15:  Magnitude  of  piecewise  constant  IBC  term, 
terms  for  a  35  yum  thick  conductor  with  cr  =  5.8  x  10^  fl 
At  =  (dash),  At  =  (dash-dot). 


^^AmAt  _  .^g  expansion 

-^m-i;  At  =  ^  (solid). 


2c 


2c 
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Figure  4-16:  Piecewise  constant  error  vs.  expansion  terms  for  various  conductor 
thicknesses  with  A  =  .005  m  and  cr  =  5.8xl0’^S7“^m“^;  I  =  17.5 //m  (solid),  I  =  35  [xm 
(dash),  I  =  70 ixm  (dash-dot). 

4.7  Verification  of  IBC  in  Two  Dimensions 

Measurement  of  the  quality  factor  of  resonant  structures  is  used  to  verify  the  multi  - 
dimensional  IBC.  This  method  is  desirable  for  two  reasons.  First,  it  is  the  very 
small  but  finite  loss  that  a  good  conductor  exhibits  that  causes  errors  in  FDTD 
simulations  where  good  conductors  are  represented  by  perfect  conductors.  Second, 
the  measurement  is  based  on  many  reflections  from  the  cavity’s  surfaces  and  should 
give  a  good  indication  of  the  accuracy  and  thus  the  benefit  of  this  type  of  IBC. 

4.7.1  Calculation  of  Quality  Factor 

An  analytical  approximation  for  the  quality  factor  of  a  resonant  cavity  is  used  for  the 
comparison  against  the  FDTD  IBC.  Using  the  nomenclature  and  derivation  found  in 
[138],  the  essential  ideas  follow. 

The  quality  factor  of  Q  is  a  measure  of  the  loss  of  a  resonant  cavity.  Let  (t) 
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represent  the  total  electromagnetic  energy  of  a  mode  n  within  a  resonant  cavity.  The 
energy  can  be  approximated  by 


WT„{t)  =  wr„(0)e  On 


(4.39) 


where  Wn  is  the  resonant  frequency  in  radians  of  the  nth  mode.  So  using  the  definition 
in  [138],  the  quality  factor  is  the  average  number  of  radians  it  takes  for  the  total 
electromagnetic  energy  to  decay  to  ^  of  its  original  value.  The  higher  the  Q  the  lower 
the  loss.  The  quality  factor  can  be  determine  from  the  average  power  dissipated  within 
the  structure  and  the  average  electromagnetic  energy  stored  within  the  structure  given 
by 


Qn  — 


Uji  X  energy  stored  in  nth  mode 
average  power  dissipated  in  nth  mode 


{Pn)  ■ 


(4.40) 


The  power  lost  in  the  resonant  structures  used  in  this  chapter  are  a  result  of  the 
current  flowing  through  the  walls  of  the  cavity.  So  to  calculate  the  Q  of  the  resonant 
structure  under  test  we  must  calculate  the  following  two  integrals: 

WT„  =  ^J^eo\En\^dv  (4.41) 

and 

Pn  =  \f^En-J:dv.  (4.42) 

As  in  [138],  since  we’re  dealing  with  very  large  conductivities  perturbation  methods 
can  be  used.  The  current  can  be  approximated  by  the  current  generated  in  a  perfect 
conductor.  So  the  power  calculation  of  (4.42)  becomes 

Pn-ll  — (4.43) 
2  Jv  a 

where  J°  represents  the  current  in  a  resonator  with  PEC  walls.  Again  using  the 
good  conductor  approximation,  most  of  the  current  flows  within  one  skin  depth  of  the 
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Figure  4-17:  Test  structure  for  IBC  validation. 


conductor,  where  the  skin  depth,  6  is  defined  by 


6  = 


(4.44) 


and  the  current  can  be  approximated  by  assuming  that  the  current  is  equivalent  to 
the  surface  current  of  a  perfect  conductor  uniformly  distributed  over  a  single  skin 
depth,  so  equation  4.43  becomes 


where  A  is  the  inside  surface  of  the  resonator  and  is  the  surface  current  due  to  the 
nth  mode  in  a  PEC  resonator.  Figure  4-17  is  a  schematic  of  the  test  structure  used 
to  validate  the  two  dimensional  IBC.  The  test  structure  is  a  rectangular  resonating 
cavity  made  with  a  good  conductor  with  conductivity,  a,  with  dimensions  (in  meters) 
a  in  the  x  direction  and  c  in  the  2  direction.  The  thickness  of  the  thin  conducting 
sheet  is  I  m  and  the  skin  depth  is  6  m.  (The  figure  is  not  drawn  to  scale  I  ^  a,  c.) 
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The  cavity  is  excited  with  TE  modes  by  using  a  line  source  located  at  {is,  kg)  on  a 
discretized  FDTD  grid  and  the  response  is  measured  with  a  software  probe  at  {ip,  kp) 
that  records  Ey{ip,  kp)  at  each  time  step. 

For  the  power  and  energy  calculations  we  assume  infinite  extent  in  the  y  direction 
and  therefore  the  energy  and  power  units  are  all  described  per  meter. 


With  Figure  4-17  in  mind  the  power  dissipated  per  meter  is  given  by 

Pn  =  2^^  1^  X  Hl{z  =  0)pd3; 

+  2^J^^\xxHl{x  =  D)\^dz.  (4.46) 

Given  a  TE  excitation  {Ey),  the  mode  can  be  described  by 

Ey{x,  z)  =  Eo  sin  kxX  sin  k^z,  (4.47) 


where  k^  =  '^,  and  k^  =  ^.  The  magnetic  fields  can  be  found  from  Maxwell’s 
equations  and  (4.47).  They  are 

k 

Hx{x,  z)  =  i — —Eo  sin  kxX  cos  kzZ  (4.48) 

UJ/lo 


(4.49) 


(4.50) 
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The  total  electromagnetic  energy  is  calculated  from  (4.41)  and  (4.47),  resulting  in 


I 


(4.51) 


Substitution  of  (4.50)  and  (4.51)  into  (4.40)  yields 

47r^[m^a^  +  p^c^]  ’ 


(4.52) 


substituting  ul^  peppLo  =  p  =  A;2  +  A;2  ^  7r2  ^  j  arrive  at 


_  {a?rn?  +  c^p^)ac 
25[m2a^  +  p'^c?] ' 


(4.53) 


The  quality  factor  of  a  good  conductor  resonator  is  a  function  of  the  mode  number, 
skin  depth,  and  the  physical  dimensions  of  the  resonator. 


4.7.2  Measurement  of  Quality  Factor  in  FDTD  Simulations 

There  are  two  ways  to  measure  the  quality  factor  in  a  FDTD  simulation.  The  first  is 
the  take  the  Fourier  transform  of  the  electric  field  measured  by  the  software  probe.  By 
measuring  the  half-power  bandwidth  about  the  resonator’s  resonant  frequencies,  the 
quality  factor  can  be  derived  [138],  [140].  However,  this  method  will  require  millions 
of  time  steps  for  very  high  Q  structures.  A  more  practical  method  relies  on  exciting 
a  single  mode  and  measuring  the  decay  in  the  mode’s  field  strength  measured  by  the 
software  probe.  In  the  lossy  resonator,  the  resonant  frequencies  are  better  described  in 
the  s  domain,  Sm,p  —  —  Since  the  energy  is  proportional  to  the  magnitude 

of  the  electric  field  squared,  equation  (4.39)  can  be  rewritten  as 


=  0)e 


(4.54) 
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MSiemans 

1/6 

Qfdtd 

Q 

%  error 

.058 

.92 

801.5 

922.0 

-13.1 

.58 

2.92 

2950.4 

2915.5 

+1.2 

5.8 

9.22 

9220.3 

9219.7 

0.0 

58.0 

29.16 

30365.7 

29155.2 

-K3.8 

580 

92.20 

100939.0 

92196.8 

-t-9.4 

Table  4.1:  Comparison  of  Qfdtd  and  analytic  Q  at  3.03  GHz 


and  the  quality  factor  becomes 


Q 


m,p 


O/y 


(4.55) 


If  a  single  mode  is  excited  in  the  FDTD  simulation,  oiTn,p  caii  be  calculate  from  two 
peaks  in  the  stored  electric  field,  given  by 


(n2  — ni)At 


(4.56) 


4.7.3  2D  Simulations  and  Results 


The  impedance  boundary  condition  is  extended  to  two-dimensions.  Using  the  test 
structure  in  Figure  4-17  the  dimensions  of  the  structure  are  changed  to  vary  the  fre¬ 
quency.  The  current  source  is  placed  at  the  center  of  the  structure  and  the  probe  is  off 
center.  A  modulated  gaussian  pulse,  modulated  at  the  TEi^i  resonant  frequency  with 
/?  =  4096,  is  used  as  the  temporal  excitation.  The  software  probe  stores  the  electric 
field  at  each  time  step.  Equation  (4.56)  is  used  to  measure  and  equation  (4.55) 
is  used  to  measure  Qfdtd-  The  simulated  results  are  compared  to  the  analytically 
derived  values  of  Q.  Tables  4.1-  4.4  show  the  results.  Each  table  represents  a  different 
frequency  and  each  table  spans  four  orders  of  magnitude  of  conductivity.  The  tabular 
results  are  shown  graphically  in  Figure  4-18. 
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MSiemans 

ijb 

Qfdtd 

Q 

%  error 

.058 

1.30 

693.5 

652.0 

+6.4 

.58 

4.12 

2019.6 

2061.6 

-2.0 

5.8 

13.04 

6521.0 

6519.3 

0.0 

58.0 

+1.2 

580 

130.39 

63140.3 

65193.0 

-3.1 

Table  4.2:  Comparison  of  Qfdtd  and  analytic  Q  at  6.06  GHz 


MSiemans 

116 

Qfdtd 

%  error 

.058 

455.4 

+13.9 

.58 

5.18 

1399.0 

1440.2 

5.8 

16.40 

4574.4 

4554.2 

58.0 

51.85 

14502.5 

14401.7 

■SsU 

580 

163.95 

45889.3 

45542.1 

+0.8 

Table  4.3:  Comparison  of  Qfdtd  and  analytic  Q  at  9.6  GHz 


MSiemans 

Ijb 

Qfdtd 

Q 

%  error 

.058 

1.84 

501.0 

461.0 

+8.7 

.58 

5.83 

1401.8 

1457.8 

-3.8 

5.8 

18.44 

4642.0 

4609.8 

+0.7 

58.0 

58.31 

14430.6 

14577.6 

-1.0 

580 

184.39 

44324.1 

46098.4 

-3.8 

Table  4.4:  Comparison  of  Qfdtd  and  analytic  Q  at  12.1  GHz 
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Figure  4-18:  2D  resonator  quality  factor  percent  error  versus  conductivity  exponent, 
X  in  5.8  X  10^f2-im-i;  3.03^77^  (solid),  6.06  GHz  (dash),  9.60 GHz  (dash-dot), 
12.10  GHz  (dash-dot-dot). 

4.8  Discussion 

The  conductivity  range  studied  in  this  chapter  spanned  four  orders  of  magnitude  from 
10“*  to  (Some  work  at  10®  was  done  in  Section  4.6.2.)  Even  the  lowest 

conductivity  used  would  be  considered  a  good  conductor  at  the  frequencies  tested 
as  defined  by  1  <<  The  conductivities,  in  of  the  most  popular  signal 

conductors,  silver,  copper,  gold  and  aluminum  are  6.14  x  10^,  5.80  x  10^,  4.10  x  10^ 
and  3.54  x  10^  and  the  conductivity  of  steel  is  3.54  x  10^  which  all  easily  fall  within 
the  range  of  conductivities  tested. 

When  a  piecewise  constant  approximation  is  used  for  the  magnetic  field,  the 
amount  of  computer  resources  (as  measured  by  number  of  expansion  terms)  is  di¬ 
rectly  proportional  to  the  conductivity.  Lower  conductivities  result  in  larger  pole 
magnitudes  and  lesser  sensitivity  to  increase  number  of  terms.  The  piecewise  con¬ 
stant  approximation  produces  an  error  that  cannot  be  overcome  by  more  expansion 
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terms.  In  fact,  the  errors  from  the  impedance  approximation  and  the  magnetic  field 
approximation  combine  to  produce  an  optimum  number  of  expansion  terms. 

With  a  piecewise  linear  approximation,  the  error  due  to  the  magnetic  field  is 
greatly  reduced  and  an  optimum  number  of  terms  is  not  evident.  The  more  terms 
used  the  more  accurate  the  IBC;  however,  the  error  levels  off  at  200  terms  in  the  test 
of  the  copper  sheet  in  Section  4.6.3. 

When  extended  to  two  dimensions,  except  for  the  lowest  conductivity  of  cr  = 
5.8  X  10'^  the  impedance  boundary  condition  works  incredibly  well  over  the 

wide  range  of  skin  depths  studied,  3  ^  185.  The  IBC  is  an  excellent  model  for  thin 
sheets  of  good  conductors  with  finite  conductivity.  The  IBC  offers  the  advantage  of 
working  over  a  wide  range  of  frequencies  with  very  little  computational  overhead.  For 
typical  conductors,  like  copper,  125  terms  per  tangential  electric  field  component  are 
needed  with  errors  in  Q  less  than  5%. 


Chapter  5 


Phase  Unwrapping  of  Synthetic 
Aperture  Radar  (SAR) 
Interferometry 

5.1  Introduction 

Airborne  and  spaceborne  Synthetic  Aperture  Radar  (SAR)  platforms  have  been  used 
for  many  years  to  study  the  earth’s  surface  [144]-[159].  When  two  radars  on  a  sin¬ 
gle  platform  or  two  passes  of  a  single  radar  map  the  same  area,  an  interferogram 
can  be  produced  from  the  difference  in  phase  measured  by  each  radar  or  pass.  An 
interferogram  is  a  pictorial  representation  of  the  phase  differences  measured  at  each 
pixel. 

Since  the  measured  phase  differences  lie  between  — tt  and  tt,  the  phase  is  said  to 
be  wrapped.  A  SAR  interferogram  contains  fringes.  These  fringes  are  the  locations 
on  the  interferogram  where  a  27r  discontinuity  exists.  The  interferogram  resembles  a 
topographical  contour  map  where  a  line  of  constant  elevation  corresponds  to  a  fringe. 
When  no  noise  is  present,  the  fringes  can  easily  be  located  and  the  data  adjusted 
by  adding  multiples  of  27r  to  produce  an  unwrapped  phase  image.  However,  real- 
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world  data  are  always  contaminated  with  noise,  hence  there  is  a  need  to  develop 
sophisticated  phase  unwrapping  algorithms. 

Successful  phase  unwrapping  is  the  key  to  the  extraction  of  DEM  (digital  elevation 
model)  from  an  interferometric  SAR  phase  image.  There  are  local  and  global  methods 
to  unwrap  phase  data.  Local  methods  unwrap  a  pixel  based  only  on  its  nearest 
neighboring  pixels;  whereas  global  methods  consider  the  whole  set  of  pixels  to  be 
unwrapped.  Phase  unwrapping  has  been  a  topic  of  much  research  [160]  -[175]. 

Local  or  global  aside,  there  are  really  two  basic  approaches  to  unwrapping  phase 
data.  The  first  is  based  on  finding  an  unwrapped  solution  such  that  the  solution’s 
first-order  partial  derivatives  in  the  x  and  y  directions  match  (or  closely  match)  the 
wrapped  first-order  partial  derivatives  or  gradients  of  the  phase  data.  Typically  noise 
is  handled  by  unwrapping  the  best  data  first  in  local  schemes,  or  in  global  schemes, 
like  the  least  squares  method,  the  data  is  weighted  to  favor  the  best  data.  The  second 
approach  moves  along  the  data  and  adds  or  subtracts  27r  when  a  fringe  is  crossed.  The 
fringes  are  found  with  a  fringe  detection  scheme  [155].  With  this  method,  one  way 
to  handle  the  noise  is  to  form  branch  cuts.  This  way  the  image  is  searched  for  phase 
inconsistencies,  in  the  form  of  residues.  The  residues  are  connected  to  form  branch 
cuts  and  phase  is  unwrapped  by  adjusting  the  integration  path  or  by  modifying  the 
fringe  information. 

There  are  local  and  global  techniques  to  phase  unwrapping  algorithms  based  on  the 
gradients  of  measured  data.  The  local  algorithms  usually  locate  one  or  more  areas  on 
the  image  that  are  considered  good  data.  To  find  these  points  one  may  use  noise  floor 
data  and  signal-to-noise  ratio  [168]  or  use  coherence  data  [178],  [166].  Unwrapping 
usually  begins  with  these  good  areas.  The  algorithms  move  from  pixel  to  pixel  and 
add  or  subtract  27r  based  on  some  criteria.  One  method  follows  the  least-gradient 
path  with  the  assumption  that  the  smallest  gradient  points  to  the  best  data  [178], 
[180]  and  the  adjustments  are  made  to  match  the  solution  gradient  to  the  wrapped 
measured  gradient.  Other  methods  use  more  neighbors  to  decide  the  value  to  add  to 
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the  unwrapped  pixel  [161],  [169].  Here  all  the  gradients  of  the  adjacent  unwrapped 
pixels  are  examined  and  the  phase  of  the  unwrapped  pixel  is  predicted  based  on  these 
gradients.  Then  an  integer  number  of  2w  is  added  to  bring  the  unwrapped  pixel  closest 
to  the  predicted  phase.  In  this  way  it  is  possible  to  have  phase  differences  larger  than 
TT,  so  this  method  can  accommodate  real  discontinuities  due  to  terrain  features.  All  of 
these  local  approaches  have  the  added  complication  of  growing  separate  unwrapped 
regions  that  must  be  joined  to  produce  the  final  product.  Since  the  best  pixels  are 
unwrapped  first,  the  likelihood  of  errors  propagating  through  the  image  is  reduced. 

Global  gradient  algorithms  are  based  on  least  squares  and  weighted  least  squares 
methods  [170].  In  this  way,  the  technique  attempts  to  find  a  solution  that  minimizes 
the  differences  between  all  of  the  solution’s  gradients  and  the  wrapped  data’s  gradi¬ 
ents.  The  least  squares  approach  is  very  desirable  because  of  the  speed  at  which  all 
the  data  can  be  unwrapped,  but,  the  method  cannot  treat  noisy  data  very  well  [170]. 
The  weighted  least  squares  allows  the  user  to  favor  the  good  data  by  applying  a  set 
of  weights  to  the  data  based  on  some  knowledge  of  the  data’s  noise  content.  These 
methods  have  been  applied  to  a  single  SAR  interferogram  in  [160];  however,  no  details 
of  the  weighting  criteria  was  presented.  Others  have  used  a  weighted  least  squares 
method  to  unwrap  a  simulated  interferogram  with  a  shear  [162],  where  the  shear  dis¬ 
continuity  was  masked  out  and  a  multi-grid  iteration  scheme  was  used  to  solve  the 
weighted  partial  differential  equation.  A  weighted  least  squares  technique  has  been 
used  unwrap  speckle  interferograms  using  a  weighting  scheme  based  on  masking  out 
the  phase  inconsistencies  or  residues  [171]. 

The  branch  cut  method  also  has  local  and  global  approaches.  So  far  only  local 
approaches  have  been  reported  in  SAR  interferometry.  The  branch  cut  method  and 
SAR  interferometry  phase  unwrapping  was  first  reported  in  [173].  Applied  to  the 
interferogram  derived  from  two  passes  of  Seasat,  the  residue  connection  was  based  on 
a  nearest  nezp'/ifeor  approach.  This  technique  works  well  with  a  low  density  of  residues; 
however,  it  breaks  down  with  a  high  density  of  residues  [173].  Once  two  residues  are 
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connected  they  are  removed  and  no  longer  considered  for  connection  to  any  other 
residues.  This  method  is  very  likely  to  leave  an  uncompensated  residue.  The  phase 
is  then  unwrapped  by  integrating  along  a  path  that  never  crosses  a  branch  cut.  This 
approach  was  briefly  reported  in  [160]  as  part  of  an  overview  of  SAR  phase  unwrapping 
techniques.  It  was  report  that  the  main  disadvantage  to  the  branch  cut  method  is 
the  propagation  of  global  errors  from  uncompensated  residues.  A  modification  to  the 
basic  branch  cut  method  in  SAR  interferometry  involves  connecting  and  removing 
residue  pairs  that  are  only  1  or  2  pixels  apart  [163].  The  remaining  residues  are 
handled  separately  and  considered  part  of  real  discontinuities  that  exist  as  a  result  of 
terrain  features. 


The  branch  cut  method  has  also  been  applied  to  speckle  interferometry,  used 
to  measure  very  small  surface  deformations  on  structures.  In  this  application,  a 
nearest  neighbor  (local)  connection  algorithm  was  used  and  reported  in  [174].  The 
first  global  branch  cut  method  is  used  to  unwrap  a  speckle  interferogram.  All  residues 
are  considered  before  making  any  branch  cuts.  In  this  way,  the  algorithm  is  based  on 
minimizing  the  total  branch-cut  length  [175]. 


Finding  efficient,  accurate  and  automatic  phase  unwrapping  algorithms  is  an  in¬ 
teresting  topic  and  is  becoming  more  important  as  more  processing  occurs  on  the  SAR 
platform.  Although  some  new  phase  unwrapping  techniques  have  been  briefly  intro¬ 
duced  recently  based  on  the  principle  of  maximum  entropy  [164]  and  multiresolution 
[165],  they  have  not  been  used  as  yet  on  SAR  data.  In  this  chapter,  the  weighted  least 
squares  approach  is  investigated  and  applied  to  SAR  interferometry.  Specifically,  the 
a  hybrid  method  that  uses  branch  cut  residues  to  weight  the  data  is  presented  and 
compared  to  weighting  schemes  based  on  coherence  data.  Also,  optimum  global  and 
local  branch  cut  algorithms  are  presented  and  applied  to  SAR  interferograms.  Data 
from  both  simulated  and  real  SAR  interferograms  are  used. 
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5.2  Least  Squares  Phase  Unwrapping  Technique 

The  least  squares  approach  attempts  to  minimize  the  difference  between  the  wrapped 
derivatives  of  the  measured  data  and  the  derivatives  of  the  solution.  It  turns  out  that 
through  a  little  manipulation  the  least  squares  problem  can  be  constructed  to  form 
a  discretized  Poisson  equation  where  the  source  term  is  calculated  from  the  wrapped 
data.  The  equation  can  be  solved  with  the  discrete  cosine  transform  (DCT)  offering 
a  fast  and  fully  automatic  algorithm. 

Consider  axi  M  x  N  matrix  of  wrapped  phase  values.  The  element  of 

the  wrapped  phase  matrix  will  be  represented  by  4>fj-  The  goal  is  to  find  the  actual 
(unwrapped)  phase,  based  on  the  measured  wrapped  phase. 

First,  a  wrapping  operator,  W  {}  ,  must  be  defined  such  that 

W{x}  =  x  +  k2Tr,  (5.1) 

where  k  is  any  integer  that  satisfies 

— TT  <  X  +  k2'K  <  TT.  (5.2) 

The  wrapped  phase  values  may  be  described  in  terms  of  the  unwrapped  phase  values 
by 

=  +  (5-3) 

where  riij  is  the  noise  associated  with  pixel  {i,j).  The  least  squares  method  uses  the 
first  order  wrapped  differences  of  the  data  given  by 

^  0<i<M-2-,  0<j<N-l 

0  otherwise 


144 


CHAPTER  5.  PHASE  UNWRAPPING  OF  SAR  INTERFEROMETRY 


and 


}  0<i<M-l-0<i<N-2) 

1  0 

otherwise  J 

(5.5) 


The  solution  (pi^j  is  the  set  of  all  such  that  the  sum  of  all  the  squares  of  the 


differences  between  the  solution  phase  differences,  4>i+\o  —  4>i,j  and  (j>ij+i  —  and 
the  measured  phase  differences,  and  Af^,  is  a  minimum.  In  other  words,  is 
the  least  squares  solution  to 


M-2N-1  M-lN-2 

z  E  - <i>ij  (5.6) 

iz=0  j=0  2r=0  j=0 

Equation  (5.6)  can  be  formulated  as  a  discretized  Poisson  equation  with  the  wrapped 
data  as  the  basis  of  the  source  term.  In  continuous  space  the  Poisson  equation  is 

y)  +  y)  =  y)  (5-7) 

where  0(a;,  y)  is  the  continuous  phase  function  and  p{x,  y)  is  a  continuous  source  term. 
The  discretize  form,  given  the  sampled  phase  values,  is 


{4>i+l,j  —  +  4>i-l,j)  +  (^ij+l  —  =  Pij  (5.8) 

where 

Pid  =  -  ^Ud)  +  (5-9) 

For  those  source  values  that  lie  on  the  edge  of  the  domain  the  phase  differences  are 
set  to  zero  as  (5.4)  and  (5.5)  indicate.  Essentially  a  Neumann  boundary  condition  is 
imposed.  Now  that  a  discrete  Poisson  equation  has  been  constructed  the  solution  can 
be  found  with  the  use  of  the  Discrete  Cosine  Transform  (DCT).  Using  cosine  series 
expansions  for  the  unknowns  leads  to  a  very  simple  solution  technique.  The  DCT 
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and  Inverse  DCT  pair,  Xij  Cm,ni  defined  as 

M-lN-l 
3 


M— liV— i  r  1  r  / 

f  E  E  cos  ^m(2i  +  l)  cos  ^n(2j  + 1)  1 

i-0  j=0  I 

Crr,.„.  =  <  0<m<M-l;  0<n<iV-l; 

0  otherwise 


(510) 


and 


^i,j  —  \ 


f  E  E  4wi{m)w2{n)Cm,n  cos  yj^m{2i  +  1)]  cos  [^n{2j  +  1)]  I 

*  m=0  n=0  ' 

0<i<M-l;  0<j<N-l', 

0  otherwise 


(511) 


where 


Wi  =  1/2,  m  =  0, 

lUi  =  1,  1  <  m  <  M  —  1, 

W2  =  1/2,  n  =  0, 

^2  =  1,  l<n<N-l, 

Substituting  the  cosine  expansions  into  (5.8)  leads  to 


(512) 


_ Phj _ 

2  (cos  +  cos 


(5.13) 


where  and  pij  are  the  discrete  cosine  transforms  of  4>ij  and  Pij. 

The  solution  may  be  found  by  performing  a  2D  DCT  on  the  source  term,  modifying 
the  transformed  source  term  by  using  the  right  hand  side  of  (5.13),  and  performing 
an  inverse  DCT  on  the  modified  source  term  to  arrive  at  the  solution  Very  fast 
DCT  algorithms  have  been  developed  based  on  the  Fast  Fourier  Transform  [143]  and 
can  be  used  to  unwrap  phase  very  quickly. 


The  disadvantage  to  the  least  squares  method  is  that  it  does  not  consider  the 
integrity  of  the  data.  If  only  a  small  localized  portion  of  the  interferogram  is  corrupted 
by  noise,  the  entire  unwrapped  image  is  effected  [170].  In  other  words,  a  local  error 
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will  propagate  throughout  the  image. 


5.3  Weighted  Least  Squares  Phase  Unwrapping  Tech¬ 
nique 

The  least  squares  solution  to  phase  unwrapping  assumes  that  the  noise  associated 
with  each  pixel  is  Gaussian  and  independent  [143].  The  assumptions  will  not  always 
apply  to  the  noise  in  a  SAR  interferogram.  Applying  a  weighting  scheme  attempts 
to  incorporate  the  knowledge  of  the  interferogram’s  pixel  integrity  when  unwrapping 
the  phase. 

5.3.1  Weighted  Least  Squares 

If  we  suspect  that  a  measured  phase  value  is  are  corrupted  with  noise  or  otherwise 
unreliable,  as  assumed  in  real  data,  we  can  weight  that  phase  value,  4'f,j  by  Wij  where 

0.0  <  Wij  <  1.0.  (5.14) 

The  best  values  are  assigned  weights  of  one  or  close  to  one  whereas  the  worst  data 
are  assigned  weights  near  zero.  Although  the  weighted  least  squares  problem  cannot 
directly  be  solved  with  the  Discrete  Cosine  Transform,  Ghiglia  et  al.  [170]  showed 
how  an  iterative  approach  can  be  constructed  to  take  advantage  of  the  DCT.  Instead 
of  a  source  term  based  solely  on  the  measured  data,  a  new  source  term  built  with  a 
weighted  sum  of  the  measured  data’s  first  derivatives  and  the  current  iteration’s  first 
derivative  is  used. 

The  first  order  differences  of  the  data  are  weighted  by  f^f  such  that 

fij  = 


(5.15) 
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and 

ffj  = 

The  portion  of  the  source  term  based  on  the  data,  pf  p  is  given  by 

pfj  =  {ftAh  -  fi-hAU,i)  +  UlAh  -  fh-iKi-i) 

and  the  contribution  from  the  current  iteration’s  solution,  \  is  given  by 


(5.16) 


(5.17) 


so  that  the  A:th  solution  is  found  by  solving  the  discretized  Poisson’s  equation  with  a 
new  source  term. 


Ph 


d  I  ^(A;- 

Pip  +  Pip 


1) 


(5.19) 


Solving  the  weighted  least  squares  method  in  this  way  is  known  as  the  Picard  Iteration 
Method  and  is  described  in  Section  5.3.3.  The  iterations  are  performed  until  a  user 
specified  number  of  iterations  are  complete.  The  method  is  still  very  fast  since  it  uses 
the  DCT  and  the  total  time  depends  on  the  number  of  iterations  performed.  Note 
that  if  all  Wij  =  1  then  the  weighted  least  squares  method  becomes  the  least  squares 
method  of  Section  5.2. 


5.3.2  Determining  Weights 

When  using  the  weighted  least  squares  method  for  unwrapping  interferograms,  some 
scheme  must  be  used  to  assign  the  weights  to  each  pixel.  In  this  section,  weighting 
schemes  based  on  coherence  data  are  described. 
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Coherence  Weighting 

The  coherence  of  the  image  is  directly  related  to  the  integrity  of  the  data.  Given  the 
complex  returns  from  both  sensors  due  to  the  same  place  on  the  earth,  and 

S2{hj),  the  coherence,  p,  of  pixel  {i,j)  is  defined  by 


V  j)S2{i,  j)S^{i,  j) 


(5.20) 


The  maximum  value  for  the  coherence  is  1  and  the  minimum  value  is  0.  Low  coherence 
can  be  the  result  of  many  factors  such  as  thermal  noise  and  processing  phase  errors. 
In  repeat-pass  spaceborne  systems,  decorrelation  can  come  from  different  atmospheric 
conditions  during  each  pass,  changes  in  the  terrain  between  the  passes,  and  from  the 
different  viewing  positions  of  the  passes[155].  Also  the  length  of  the  baseline  is  a 
critical  factor.  If  the  baseline  is  too  large  there  is  a  complete  loss  of  coherence  and 
the  corresponding  phase  data  is  useless  for  terrain  height  inversion. 

The  coherence  data  is  a  logical  choice  for  a  weight.  When  coherence  weighting  is 
used,  Wij  =  Cij  where  Wij  is  the  weight  assigned  to  pixel  {i,j)  and  Ctj  is  the  coherence 
of  pixel  {i,j)  of  equation  (5.20). 


Modified  Coherence  Weighting 

A  variant  of  the  straight  coherence  weighting  scheme  is  to  set  a  coherence  threshold, 
Cthr,  in  order  to  give  the  maximum  weight  to  good  data.  Table  5.1  shows  the  algorithm 
used  to  produced  modified  coherence  weights.  Setting  Cthr  =  0  is  equivalent  to  no 
weighting  {i.e.  a  straight  least  squares  method)  and  setting  Cthr  =  1  is  equivalent  to 
straight  coherence  weighting. 


5.3.3  Picard  Iteration  method 

As  stated,  the  Picard  iteration  method  is  used  to  solve  the  weighted  least  squares 
problem.  Following  the  description  found  in  [170],  starting  with  the  over  determined 
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Set  threshold  Cthr 
If  Cij  >  Cthr  then 

Wij  =  1 

else 

endif 


Table  5.1:  Algorithm  for  weight  assignments  based  on  coherence  data. 


problem  given  by  the  matrix  equation 


Ax  = 


b, 


(5.21) 


the  least  squares  solution  is  found  by  solving 

A  Ax  =  A  b, 


(5.22) 


where  x  is  the  solution  vector  containing  MN  phase  values  and  6  is  a  vector  containing 
N{M  —  1)  +  M{N  —  1)  wrapped  phase  differences.  Let  P  =  A  A  (the  discrete 

=T_ 

Laplacian),  p  =  A  b  (the  source  term),  and  (/>  be  the  phase  vector;  then  the  least 
squares  solution  of  the  phase  unwrapping  problem  is  given  by 


P^  =  p. 


(5.23) 


However,  if  we  apply  a  set  of  weights  to  the  wrapped  phase  differences  b  then  (5.21) 
becomes 


WAx  =  Wb, 


(5.24) 


and  the  solution  is  found  by  solving 


A^w'^W  Ax  =  AW^W  b. 


(5.25) 
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Defining  Q  =  W  A,  c  =  W  b,  and  4>  as  the  phase  vector  again,  the 

system  of  equations  to  solve  is 

S  =  c.  (5.26) 

The  discrete  cosine  transform  cannot  be  used  to  solve  (5.26)  as  written.  The  Picard 
method  rewrites  Q  as  the  sum  of  the  discrete  Laplacian,  P,  and  a  difference  matrix, 
D  such  that 

Q  =  Tt^.  (5.27) 

Substitution  of  (5.27)  into  (5.26)  and  a  little  algebra,  leads  to 


P(j)  =  c  —  D(j). 

Equation  (5.28)  can  be  solve  iteratively  using  the  DCT  and 

P^k+i  =  c-  D^k, 

since  (5.29)  is  the  discrete  Poisson  equation, 

P4*k+1  ~  Pki 

where 


from  equation  (5.19). 


Pk 


c  —  D((>k 

/ .  + 

I  Hi,j 


(5.28) 


(5.29) 


(5.30) 


(5.31) 


A  disadvantage  of  the  Picard  iteration  method  is  that  convergence  is  not  guaran¬ 
teed.  The  user  selects  the  number  of  iterations  and  the  algorithm  terminates  when 
that  iteration  number  is  reached. 
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5.4  Branch  Cut  Method 

The  branch  cut  method  is  very  different  from  the  least  squares  approach.  Instead 
of  forming  a  set  of  equations  to  solve,  the  solution  is  derived  from  integrating  the 
fringes.  This  approach  is  intuitively  appealing,  but  offers  its  own  set  of  difficulties. 

Unwrapping  phase  using  the  branch  cut  method  begins  by  constructing  a  fringe 
map  from  the  wrapped  data.  A  fringe  map  records  the  locations  of  the  2%  discontinu¬ 
ities  call  fringes.  The  wrapped  data  is  unwrapped  by  moving  through  the  image  and 
integrating  the  fringes.  Whenever  a  fringe  occurs,  2ir  is  either  added  or  subtracted 
to  the  data  depending  on  the  value  of  the  fringe.  In  the  presence  of  noise,  erroneous 
fringes  occur  and  cause  errors  to  propagate  through  the  image.  If  phase  inconsisten¬ 
cies  can  be  found  then  steps  can  be  taken  to  reduce  the  effects  of  the  noise.  Since  we 
know  that  there  should  be  an  equal  number  of  positive  and  negative  fringes  if  we  travel 
around  a  closed  contour,  we  can  search  for  phase  inconsistencies  by  integrating  the 
along  a  closed  contour.  If  the  integration  is  non-zero  then  the  contour  must  enclose 
residues.  The  presence  of  residues  means  the  unwrapped  solution  is  path  dependent. 
To  isolate  the  residues,  we  integrate  along  the  smallest  contour  possible  i.e.  between 
four  adjacent  data  points.  In  order  to  remove  the  ambiguity  of  the  integration  path 
the  residues  of  opposite  sign  are  connected  to  form  branch  cuts.  The  phase  fringe 
map  is  then  modified  by  incorporating  branch  cuts  and  the  phase  is  unwrapped  in  the 
normal  fashion.  There  is  no  global  error  with  proper  treatment  of  the  single  residues. 
The  disadvantage,  however,  is  that  manual  operation  is  usually  required  to  complete 
the  branch  point  connections. 

We  begin  with  another  operator  dij  {}  defined  as 

d,j  =  {dlj,cl?,jf ,  (5.32) 

'C  -  C-i., 

27r 


where 


{c} 


(5.33) 
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and 

{C}  =  [  ■  (5.34) 

The  square  brackets,  [],  represent  rounding  to  the  nearest  integer.  The  dij  {}  operator 
will  find  any  27r  discontinuities  as  one  moves  across  the  wrapped  data  in  the  vertical 
(along  rows,  dfj)  and  horizontal  (along  columns,  d^j)  directions. 

If  no  noise  is  present  and  one  phase  pixel  value  is  known,  then  the  phase  at  all 
other  pixels  can  be  determine  by  integrating  the  number  of  discontinuities  along  any 
path  between  the  known  pixel  and  the  pixel  of  interest.  To  illustrate,  let  be  the 
known  phase  value  located  at  {io,jo)i  then 

^i,j  ~  4^io,jo  A  k2fTT  (5.35) 

where 

k  =  Yj^id-ds  (5.36) 

and  s  is  any  path  from  (ig,  jo)  to  (i,  j).  If  s  is  a  closed  path  then  k  should  equal  zero  i.e. 
we  should  end  where  we  begin.  Otherwise,  there  must  be  noise  sources  located  inside 
the  closed  contour  s.  With  these  noise  sources  or  residues,  the  simple  integration 
technique  of  (5.35,  5.36)  will  propagate  errors  throughout  the  image.  Furthermore, 
the  results  will  depend  on  the  path  taken.  Clearly  the  noise  must  be  considered  such 
that  a  unique  solution  can  be  found. 

Recall  that  the  branch  cut  method  starts  by  locating  the  27r  phase  discontinuities 
(or  fringes)  with  the  d{}  operator.  A  fringe  map  is  constructed  that  marks  the  loca¬ 
tion  of  the  fringes.  Next,  the  residues  are  found  by  integrating  the  phase  differences 
of  four  adjacent  phase  data  points  using 

resij  =  -h  ~  (5.37) 

If  resij  =  -1  we  have  a  positive  residue  and  if  resij  =  1  we  have  a  negative  residue. 
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These  definitions  are  consistent  with  the  Cauchy  Theorem  from  complex  mathematics. 
After  all  the  residues  are  located,  branch  cuts  are  drawn  from  the  positive  residues  to 
the  negative  residues.  The  branch  cut  represents  a  2tt  discontinuity.  The  direction  of 
the  discontinuity  depends  on  the  direction  of  the  branch  cut.  For  example  if  a  branch 
cut  is  pointing  down  then  the  branch  cut  represents  a  +27r  jump  when  integrating 
across  the  branch  cut  from  left  to  right.  Finally  the  fringe  map  is  updated  to  reflect 
the  branch  cut  information  and  the  solution  that  follows  is  now  path  independent. 

Although  the  concept  is  straightforward,  the  difficulty  of  the  branch  cut  method 
becomes  apparent  when  connecting  branch  points.  Since  the  residues  must  be  con¬ 
nected  in  pairs  {i.e.  positive  to  negative)  and  there  is  no  guarantee  that  the  interfer- 
ogram  will  contain  an  even  number  of  positive  and  negative  residues,  it  is  very  likely 
that  single  poles  will  remain  after  a  branch  cut  algorithm  has  been  run.  Incorrect 
treatment  of  these  single  poles  will  result  in  errors. 

The  question  is,  given  a  set  of  positive  and  negative  residues,  what  is  the  best 
way  to  connect  them?  Since  noise  will  create  pairs  of  residues  (one  positive  and 
one  negative)  the  connection  method  must  favor  connections  of  residues  close  to  each 
other.  If  all  the  residues  reside  completely  within  the  edges  of  the  SAR  interferogram, 
the  connections  can  be  made  in  a  fairly  straight  forward  manner  and  any  errors  that 
result  tend  to  be  local  in  nature.  However,  in  a  SAR  interferogram  the  data  is 
truncated  and  some  of  the  residues’s  corresponding  matches  are  not  included  in  the 
SAR  interferometry  data  set.  How  a  connection  scheme  deals  with  this  problem  will 
determine  the  algorithm’s  accuracy. 

So  far  SAR  interferometry  phase  unwrapping  has  been  limited  to  the  local  nearest 
neighbor  connection  approach  [173], [160].  This  method  can  lead  to  uncompensated 
residues  and  has  great  difficulty  when  the  residue  density  is  high.  To  alleviate  some 
of  these  problems  a  global  approach  can  be  used  that  considers  all  residues  before  the 
connections  are  made;  in  other  words  optimize  the  branch  cut  connections.  In  this 
case,  optimize  means  to  find  the  set  of  branch  cuts  that  has  the  shortest  total  branch 
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cut  length. 

5.4.1  Minimum  Branch  Cut  Algorithm 

Given  a  set  of  positive  and  negative  residue  locations  from  an  interferogram,  the 
minimum  branch  cut  algorithm  optimizes  the  set  of  connections  by  finding  the  set 
that  minimizes  the  total  branch  cut  length.  With  a  set  of  N  residue  pairs,  the 
algorithm  minimizes  the  quantity  defined  by 

&  =  E  f .  (5.38) 

k=l 

The  term,  l^,  is  the  length  of  the  kth  branch  cut  given  by 

=  \/iik-iiy  +  ijk-fk)^,  (5.39) 

where  i'^kOk)  pixel  locations  of  the  Hh  positive  and  negative 

residues  of  branch  cut  k. 

The  optimization  algorithm  is  essentially  the  same  as  the  transportation  problem 
found  in  Graph  Theory  or  linear  programming.  In  the  transportation  problem,  there 
are  m  locations  demanding  a  commodity  and  there  are  n  different  transportation 
means  to  satisfy  all  of  the  demands.  The  problem  is  to  minimize  the  total  cost  given 
that  the  unit  cost  of  using  the  ffh  transportation  vehicle  to  ship  the  commodity  to 
the  jth  location  is  Cjj.  Appendix  C  gives  an  overview  of  the  transport  problem  and 
how  it  can  be  used  to  find  the  optimum  set  of  branch  cuts. 

Before  the  transport  method  can  be  used,  the  set  of  residues  must  be  prepared. 
Given  an  interferogram,  the  first  step  is  to  find  all  the  residues  which  we’ll  call  true 
residues.  Next,  image  (image  meaning  opposite  in  sign)  residues  are  placed  at  the 
boundary  for  every  true  residue  within  a  specified  distance  from  the  edge  of  the 
interferogram.  This  distance  should  be  larger  than  the  expected  mean  branch  cut 
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length  [175].  In  this  work  about  twice  the  mean  length  is  used.  Once  all  of  the 
image  residues  are  placed,  the  positive  and  negative  residues  are  counted.  If  the 
count  is  unequal  additional  residues  are  added  to  the  boundary  which  we  will  call 
false  residues.  There  must  be  equal  numbers  of  positive  and  negative  residues  in 
order  to  perform  the  optimization.  It  is  important  to  note  that  the  problem  must 
be  constructed  so  that  false  residues  are  never  allowed  to  pair  up  with  true  residues; 
they  can  only  be  paired  with  image  residues. 

An  initial  cost  matrix  is  built  with  the  columns  representing  the  negative  residues 
and  the  rows  representing  the  positive  residues.  The  matrix  element  (*,j)  is  filled 
with  distance  between  the  ith  positive  residue  and  the  jth.  negative  residue  using 
(5.39).  (Actually  the  square  of  length  is  used  because  the  transport  algorithm  used 
demands  integer  cost  values.) 

Now,  the  matrix  is  modified  to  take  into  account  both  the  boundary  image  residues 
and  the  false  residues.  First,  all  matrix  locations  that  represent  connections  between 
any  image  residue  with  any  other  opposite  image  residue  or  false  residue  are  set  to 
zero.  In  this  way  the  algorithm  can  reduce  certain  types  of  errors  by  being  able  to 
connect  any  boundary  residue  to  any  other  boundary  residue  of  opposite  sign.  Second, 
to  ensure  no  false  boundary  residues  are  connected  to  any  true  residue,  sufficiently 
high  cost  must  be  inserted  into  the  cost  matrix  at  the  appropriate  locations.  This 
can  easily  be  done  by  placing  the  false  residues  at  fictitious  locations  far  away  from 
any  actual  interferogram  pixel. 

For  each  column,  the  minimum  length  within  the  column  is  subtracted  from  every 
entry  in  the  column.  This  step  places  a  zero  in  the  row  of  the  column  that  represents 
the  closest  positive  residues  to  the  ^‘th  negative  residue.  Then  for  each  row  the 
minimiiTTi  length  within  the  row  is  subtracted  from  every  entry  in  the  row.  This  step 
places  a  zero  in  the  column  of  the  row  that  represents  the  closest  negative  residues 
to  the  ith  positive  residue.  Now  we  have  what  Buckland  [175]  called  the  reduced  cost 
matrix. 
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Figure  5-1:  Pole  connection  ambiguity. 
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Figure  5-2;  Effect  of  missing  one  residue  in  nearest  neighbor  search. 

At  this  point,  the  reduced  cost  matrix  is  passed  to  a  subroutine  based  on  a  mini¬ 
mum  cost  transport  problem  found  in  [181].  The  routine  will  return  the  optimum  set 
of  branch  cuts  (See  Appendix  C). 

5.4,2  Branch  Cut  Limitations 

In  general  the  residue  connections  are  made  to  keep  the  branch  cut  lengths  as  small 
as  possible.  However,  several  difficulties  occur  even  though  an  optimizing  algorithm 
is  used. 

First,  there  is  ambiguity  in  the  residue  connections.  Consider  a  set  of  four  residues 
shown  in  Figure  5-1  where  the  circles  and  squares  represent  positive  and  negative 
residues  respectively.  There  are  two  equally  appropriate  connections  based  on  the 
shortest  possible  paths.  Where  one  is  correct  and  the  other  is  wrong  and  the  wrong 
choice  produces  a  local  error.  This  localized  phase  error  may  look  like  small  square 
discontinuities  in  a  relatively  smooth  image.  There  is  nothing  that  can  be  done  to 
prevent  this  type  of  error. 

Second,  algorithms  based  on  connecting  the  nearest  neighbors  invariably  will  miss 
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Figure  5-3:  Effect  of  erroneous  residue  placed  on  boundary. 

a  residue  and  force  a  connection  between  two  residues  relatively  far  apart  as  shown 
in  Figure  5-2.  This  type  of  error  (assuming  the  right  figure  is  the  correct  connection) 
produces  discontinuities  in  the  unwrapped  phase  proportional  to  the  erroneous  branch 
cut  length.  The  minimum  cost  algorithm  essentially  eliminates  this  type  of  branch 
cut  error. 

Third,  there  are  errors  that  occur  owing  to  the  treatment  of  the  edges.  With  a 
minimum  cost  algorithm  image  residues  are  placed  on  the  edges  for  those  residues  close 
to  the  edge  to  compensate  for  the  fact  that  the  residue’s  correct  pair  may  have  been 
removed  by  the  truncation  of  the  SAR  image.  An  error  may  occur  if  two  erroneous 
image  residues  are  added  as  shown  on  the  left  side  of  Figure  5-3.  Connections  between 
image  residues  must  be  allowed  as  shown  on  the  right  side  of  Figure  5-3.  Errors  like 
these  can  cause  large  errors  depending  on  the  locations  of  the  two  erroneous  residues. 
When  there  are  areas  of  high  residue  density  that  extend  between  two  edges  of  the 
interferogram,  even  the  optimum  connection  may  produce  errors  of  this  type. 

The  last  error  that  may  occur  also  stems  from  the  placement  of  image  residues  at 
the  boundary.  Clearly  we  will  not  place  an  image  residue  on  the  boundary  for  every 
residue  in  the  SAR  interferogram,  but  some  distance  (threshold)  from  the  edge  must 
be  set.  If  the  distance  is  too  small,  errors  appear  as  shown  in  the  top  figure  in  Figure 
5-4.  But  if  the  distance  from  the  edge  is  increased  slightly  more  image  residues  are 
placed  at  the  boundary  and  the  total  branch  cut  length  is  reduced. 
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Figure  5-4:  Effect  of  too  few  image  residues  on  the  image  boundary. 


We  can  never  be  certain  that  the  set  of  branch  cut  connections  is  the  correct  set. 
However,  since  noise  produces  pairs  of  residues  (positive  and  negative),  the  shortest 
branch  cut  lengths  are  more  likely  candidates  for  the  correct  connections.  Also,  the 
erroneous  connections  usually  only  produce  small  local  errors  in  the  unwrapped  phase. 

One  serious  limitation  of  the  global  optimum  branch  cut  method  is  the  CPU 
time  necessary  to  find  the  optimal  connection  set.  Since  the  solution  time  is  directly 
proportional  to  the  number  of  residues  and  proportional  to  the  square  of  the  number 
of  connections,  it  can  take  any  where  from  seconds  to  hours  to  find  the  solution. 


5.5  Local  optimum  branch  cut  method 

To  overcome  the  potential  for  long  solution  times  a  local  optimum  technique  is  intro¬ 
duced.  After  image  residues  are  added,  a  nearest  neighbor  connection  scheme  is  used. 
Starting  at  the  center  of  the  image,  the  residue  pairs  that  are  only  one  pixel  apart  are 
connected  and  removed,  then  those  pairs  two  pixels  apart,  then  three  etc.  As  residues 
are  left,  they  are  connect  to  the  closest  residue  leftover.  No  residues  are  left  uncon¬ 
nected.  Then  a  branch  cut  length  threshold  is  set.  Starting  with  the  longest  branch 
cut  that  is  greater  that  the  desired  threshold,  a  subset  of  the  interferogram  residues 
that  includes  the  long  branch  cut  is  reconnected  using  the  optimum  algorithm.  If  a 
branch  cut  cannot  be  shortened  in  the  first  pass,  a  larger  subset  of  residues  is  used.  If 
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Set  all  Wij  =  1 
If  resij  7^  0  then 

Wij  =  0 
Wij+i  =  0 
Wj+ij  =  0 

=  0 

endif 

Table  5.2:  Algorithm  for  weight  assignment  based  on  residue  identification. 


the  cut  is  still  not  shortened,  the  algorithm  moves  to  the  next  long  branch  cut.  The 
algorithm  terminates  when  attempts  have  been  made  to  fix  all  branch  cuts  longer 
than  the  threshold  length. 


5.6  Hybrid  method  -  weighted  least  squares  with 
residue  weights 

The  residue  weighting  scheme  is  a  new  method  based  on  residues  found  in  the  branch 
cut  unwrapping  method  to  target  bad  phase  values.  The  idea  is  that  if  a  residue 
exits  then  one  or  more  of  the  four  pixels  surrounding  the  residue  are  bad,  so  let’s  not 
use  them.  Table  5.2  shows  the  algorithm  used  to  produced  weights  based  on  residue 
identification. 


5.7  Comparison  of  Algorithms  on  Simulated  Data 

In  order  to  begin  a  comparison  of  the  phase  unwrapping  algorithms,  a  simulated 
interferogram  is  generated.  By  adding  noise  to  the  simulated  interferogram,  applying 
a  phase  unwrapping  process,  and  comparing  against  the  noise-free  unwrapped  data, 
we  can  begin  to  analyze  the  different  algorithms.  Two  measurements  were  taken:  the 
mean  phase  error  and  the  phase  error  standard  deviation.  The  mean  phase  error  is 
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defined  as 

kh 

and  the  phase  error  standard  deviation  is 


(5.40) 


0-0  (AV')  = 


^  ^  A-tp)  -  ^o{i,  j)  -  (j){A'ip)]‘^ 

i=i  j=i 


MN 


(5.41) 


where  is  the  root  mean  square  (rms)  phase  noise  in  degrees,  Aip)  is  an 

unwrapped  phase  value  of  a  noisy  interferogram  and  4>o{i,j)  is  the  unwrapped  phase 
data  with  no  noise  present  (A'0  =  0)  used  as  a  reference.  M  and  N  are  the  number 
of  pixel  rows  and  columns  respectively. 


5.7.1  Generation  of  Simulated  SAR  Interferogram 

In  this  work,  a  simulated  SAR  interferogram  is  generated  using  the  SAR  configuration 
is  given  in  Table  5.3.  The  interferogram  is  constructed  by  first  generating  a  simulated 
mountain  terrain.  Then  using  the  SAR  configuration,  the  phase  of  the  signal  at  each 
sensor  is  calculated  for  each  pixel  using 

=  2/conj,  a;  =  1,  2;  (5.42) 

where  ko  is  the  free  space  wave  number  and  Vij  is  the  distance  from  the  sensor  to  the 
simulated  terrain  height  for  pixel  The  interferogram  phase  value  at  pixel  {i,j) 
is  found  by  wrapping  the  phase  difference  (l)T  —  4)]^ 

Figure  5-5  is  a  gray  scale  plot  of  the  simulated  SAR  interferogram  with  no  noise 
and  Figure  5-6  is  a  gray  scale  plot  of  the  simulated  SAR  interferogram  with  A^  =  55° 
rms  noise.  These  two  figures  represent  the  range  of  rms  noise  used:  Afip  =  0°  ^  55°. 
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B 

103.0  m 

baseline 

a 

30° 

baseline  angle 

X 

.0567  m 

radar  wavelength 

9 

30° 

looking  angle 

Ho 

400  km 

sensor  height 

Table  5.3:  SAR  parameters  for  simulated  SAR  interferogram. 


Figure  5-5:  Simulated  SAR  interferogram  without  noise. 
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Figure  5-6:  Simulated  SAR  interferogram  with  Aip  =  55°  rms  noise. 


5.7.2  Addition  of  noise 


The  noise  is  added  by 

HiJ)  =  M'i'J)  +  (5.43) 

where 

<l>N{i,j)  =  A/p7r(1.0  -  2.0r„);  (5.44) 

where  Np  is  an  input  noise  parameter  to  control  the  amount  of  noise  and  r„  is  a 
random  variable  that  is  uniformly  distributed  between  0  and  1.  Therefore, 


(5.45) 

(5.46) 


3 


(5.47) 
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AV’ 

A'lp 

Np 

Alp 

Np 

IQi 

HH 

5 

0.0481 

10 

0.0962 

■a 

KOI 

IQI 

itlbVJ 

25 

0.241 

0.289 

Msm 

lifcfcfcl 

40 

0.385 

Iggjggl 

60 

0.577 

mm 

m 

liBQggl 

Table  5.4:  Np  settings  for  rms  phase  noise. 


^4>n 

A'tp 


NpTT 

x/3’ 

Np{m) 

Vs 


(5.48) 

(5.49) 


Table  5.4  lists  the  values  of  Np  necessary  produce  various  values  of  rms  phase  noise, 
A'tp. 


5.7.3  Application  and  Analysis  of  Techniques  on  Simulated 
Data 

In  this  section  the  least  squares,  weighted  least  squares  and  branch  cut  phase  unwrap¬ 
ping  techniques  are  applied  to  the  simulated  SAR  interferogram  of  the  last  section. 
Figures  5-7,  5-8,  and  5-9  are  plots  of  the  mean  and  standard  deviation  of  the  phase  er¬ 
ror  for  the  three  unwrapping  techniques  describe  earlier  in  this  chapter.  The  weights 
used  for  Figure  5-8  were  base  on  the  residue  mask  of  Table  5.2.  Table  5.7.3  lists  the 
number  of  residues  found  in  the  simulated  SAR  interferograms  for  each  rms  noise 
level  used.  Since  the  additive  noise  does  not  produce  many  phase  inconsistencies 
(i.e.  residues),  until  the  rms  phase  noise  reaches  A'0  =  45°,  we  would  not  expect  the 
weighted  least  squares  or  the  branch  cut  method  to  improve  the  error  until  Aif;  =  45°. 
The  results  show  that  both  methods  greatly  reduce  the  phase  errors. 

Figure  5-10  shows  the  standard  deviation  for  Axp  =  35°,  45°,  55°  versus  the  it¬ 
eration  number  of  the  weighted  least  squares  technique  using  the  Picard  iteration 
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Noise  Level 

Positive 

Residues 

Negative 

Residues 

Total 

Residues 

0° 

0 

0 

0 

15° 

0 

0 

0 

25° 

0 

0 

0 

35° 

17 

17 

34 

45° 

481 

484 

966 

55° 

2486 

2497 

4983 

Table  5.5:  Number  of  residues  in  simulated  SAR  interferograms. 


method.  The  mean  error  is  essentially  zero  for  all  iterations  (and  methods)  due  to 
the  nature  of  noise  added  and  not  shown.  However,  the  phase  error  standard  devia¬ 
tion  is  greatly  reduced  as  the  iteration  number  increases.  Fifty  iterations  is  sufficient 
to  reach  the  minimum  error. 

Errors  in  phase  lead  directly  to  errors  in  inverted  height.  These  relationships, 
found  in  [146],  are 


T  /  A  ,s  pAsin^  T/  A 

in 

(5.50) 

fTi(AVi)=  ,<'*(Ai/') 

47rcos(0-a)  ^ 

(6.51) 

Using  the  parameter  data  set  in  Table  5.3,  A  =  .0567  m,  $  =  30°,  a  =  30°,  {B  = 
103  m,  and  p  =  462  km  equation  (5.51)  becomes  (7ft(A'0)  =  10.1a^{A'ip)  m,  so  to 
keep  ah  <  10  m  then  the  least  squares  method  can  only  tolerate  approximately  37° 
of  rms  phase  noise  while  both  the  weighted  least  squares  and  branch  cut  method  can 
tolerate  as  much  as  55°. 


5.8  Height  Inversion 

Figure  5-11  shows  the  nominal  configuration  of  an  interferometric  SAR  set  up  used 
to  invert  terrain  height.  The  phase  difference  between  the  two  sensors,  and  S2, 
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Least  Squares 


Figure  5-7:  Least  squares  unwrapping  phase  error  versus  phase  noise,  mean  (dia¬ 
mond),  standard  deviation  (triangle). 


Weighted  Least  Squares 


Figure  5-8:  Weighted  least  squares  unwrapping  phase  error  versus  phase  noise,  mean 
(diamond),  standard  deviation  (triangle). 
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Branch  Cut 


Figure  5-9:  Branch  cut  unwrapping  phase  error  versus  phase  noise,  mean  (diamond), 
standard  deviation  (triangle). 


Figure  5-10;  Weighted  least  squares  unwrapping  phase  error  standard  deviation  versus 
iteration  number:  Atp  =  35°,  (diamond);  A'lp  =  45°,  (triangle);  Aip  =  55°,  (star). 
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Figure  5-11:  SAR  set  up  for  terrain  height  inversion. 


can  be  determined  from  the  geometry  of  the  SAR  configuration.  If  the  separation 
between  the  two  sensors,  called  the  baseline  {B),  is  much  smaller  than  the  slant  range 
(p)  to  the  target  at  P  i.e.  {B  «  p)  ,  then  the  phase  difference,  A(j)  ,  is  determined 
by  (5.52).  If  the  phase  difference  is  known,  then  9  can  be  found  from  (5.52)  and  the 
terrain  height,  h,  determined  from  (5.53). 

A'jr 

A(f)  — —B  sm{9  —  a).  (5.52) 

A 

h  =  H-pcosd.  (5.53) 

The  phase  difference  can  be  calculated  from  the  complex  radar  images  from  the  two 
sensors,  given  by 


=  arctan5i(i,  j). 


(5.54) 
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Data  Size 

2304x2048 

Horizontal  Resolution 

19.5  meters 

Vertical  Resolution 

19.5  meters 

Baseline 

233.5  m 

Range 

820000  m 

Wavelength 

0.0566  m 

Off-nadir  angle  {9) 

24° 

a 

-6.5° 

Table  5.6:  ERS  SAR  parameters. 


5.9  Application  and  Analysis  of  Techniques  on  Real 
Data 


Although  the  results  of  the  simulated  data  show  that  the  branch  cut  method  and  the 
weighted  least  squares  offer  greatly  improved  results  over  the  least  squares  method, 
the  simulated  interferograms  only  loosely  approximate  a  real  SAR  image.  In  this 
section,  the  phase  unwrapping  techniques  are  applied  to  an  actual  SAR  interferogram. 
The  height  is  inverted  and  compared  to  the  ground  truth. 


5.9.1  ERS  -1  Data 

An  interferogram  along  with  coherence  data  from  the  European  Remote  Sensing 
(ERS  -1)  satellite  was  obtained.  The  ERS  data  contains  2048  points  in  the  range 
direction  and  2304  points  in  the  azimuthal  direction.  The  interferogram  is  an  image 
of  Phillip  Smith  Mountains  in  Alaska.  The  SAR  configuration  is  contained  in  Table 
5.6.  The  ground  truth  was  obtained  from  the  U.  S.  Geological  Survey  data  of  Alaska, 
specifically  the  data  from  Phillip  Smith  Mountains  -W.  A  description  of  the  ground 
truth  data  can  be  found  in  Table  5.9.1. 


5.9.  APPLICATION  AND  ANALYSIS  OF  TECHNIQUES  ON  REAL  DATA 


169 


Data  Size: 

1201x601 

N-S  Resolution  x: 

3  arc-seconds 

E-W  Resolution  y: 

6  arc-seconds 

Lower  Left  Corner  (SE): 

68°  N 149°  fU 

Lower  Right  Corner  (NE): 

69° N  U9°W 

Upper  Left  Corner  (SW): 

68°  N  150°W 

Upper  Right  Corner  (NW): 

69°  N  150°W 

Minimum  height: 

370  m 

Maximum  height: 

2315  m 

Horizontal  accuracy: 

130  m 

Vertical  accuracy: 

30  m 

Table  5.7:  DEM  parameters. 


5.9.2  Registration  Process 


In  order  to  verify  the  phase  unwrapping  algorithms,  the  inverted  height  from  an  ERS-1 
interferogram  is  compared  to  the  ground  truth  DEM  data  from  the  U.  S.  Geological 
Survey.  Since  the  SAR  and  DEM  images  are  offset  in  the  range  and  azimuthal 
directions,  offset  by  a  rotation  angle,  and  differ  in  resolution,  the  SAR  and  DEM 
images  must  be  registered  before  any  comparisons  can  be  made.  The  registration 
process  includes  all  steps  necessary  to  synchronize  the  two  images  so  that  a  point  by 
point  compare  can  be  done. 

Throughout  the  registration  process,  the  root  mean  square  height  error  will  be 
minimized  to  find  the  registration  parameters.  The  root  mean  square  height  error  is 
defined  as 

M  N 

■Lrms  _  X  ^ 

^^error  Z-/ 

i=l j=i 


[hoEMihj)  —  hsAnihj) 


h  12 


MN 


(5.55) 


where 


I  _  ^  [hDEM{i-,3)  ~  hsAnih  j)] 

•terror  /  ^  ^ 

i=l  j=l 


MN 


(5.56) 


and  hoEMiiij)  is  the  DEM  height  interpolated  from  the  DEM  data  at  the  {x,  y)  loca¬ 
tion  of  the  SAR  pixel  {i,j)  and  hsARihj)  is  the  inverted  height  from  the  unwrapped 
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DEM  resolution  along  azimuth 
DEM  resolution  along  range 
SAR  resolution  along  azimuth 
SAR  resolution  along  range 

i^di  jd) 
(*S)  Js) 

offsetx 

offsety 
offset  angle 


^^DEM 

I^Vdem 

A.XSAR 

A^Vsar 

DEM  reference  pixel  location 
SAR  reference  pixel  location 
jd^XoEM  —  js^XsAR 
id^VDEM  —  h^VSAR 
Gos 


Table  5.8:  Registration  definitions 


interferogram  pixel  located  at  (i,i)  after  registration. 

Table  5.8  is  a  list  of  the  parameters  used  in  the  registration  process.  Furthermore, 
the  following  direction  standards  are  used.  The  azimuthal  direction  is  the  x  direction 
and  the  range  direction  is  the  y  direction.  In  addition,  the  interferogram  data  is  in 
matrix  form  where  i  indicates  the  row  and  j  indicates  the  column  and  i  =  1,  j  =  1 
is  located  in  the  lower  left  corner  of  any  displayed  data.  This  means  moving  from 
column  j  to  j  +  1  is  the  x  direction  and  moving  from  row  i  to  i  +  1  is  the  y  direction. 

The  entire  DEM  data  set  measures  601  x  1201  {i.e.  rows  x  columns).  The  x 
direction  contains  1201  data  points  separated  by  3  arc-seconds  and  the  y  direction 
contains  601  data  points  separated  by  6  arc-seconds.  The  x  direction  goes  from 
South  to  North  and  the  y  direction  goes  from  East  to  West.  Each  row  of  1201  points 
represents  a  line  of  constant  longitude  and  each  column  of  601  points  represents  a  line 
of  constant  latitude.  As  such,  Axdem  is  constant  throughout  the  image  but  AyuEM 
is  a  function  of  the  column  number.  The  resolution  decreases  as  the  column  number 
increases;  however,  the  change  is  less  than  .005%  from  column  to  column  and  less 
than  5%  over  the  entire  image. 

To  begin  the  registration  process,  the  height  must  be  inverted  from  the  SAR  data 
and  then  forshortened.  Forshortening  adjusts  the  position  of  each  pixel  along  the 
range  direction  to  compensate  for  the  fact  that  range  pixels  are  not  necessarily  in 
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order  as  far  as  the  ground  dimension  is  concerned. 

The  next  step  is  to  pick  out  the  same  feature  in  the  DEM  and  SAR  images  and 
label  those  pixels  (idjd)  and  {isjs)  respectively.  In  this  work,  the  DEM  reference 
pixel  does  not  move  while  the  SAR  reference  pixel  is  both  changed  and  moved  in 
order  to  minimize  Subsets  of  both  data  are  used  for  this  registration.  The 

subset  DEM  actual  image  size  must  be  greater  than  the  subset  SAR  image  size  so 
that  the  SAR  image  may  be  shifted,  rotated  and  compared  without  moving  off  of  the 
DEM  image.  An  initial  guess  at  the  angle  offset,  9os,  is  also  chosen. 

With  9os,  {id,jd),  and  (isJs)  as  a  starting  point,  the  SAR  reference  pixel  is  shifted 
around  {is,js)  until  a  minimum  found.  Next,  small  changes  to  9os  are  made 

until  a  minimum  is  found  again.  Since  the  SAR  data  has  a  finer  resolution 

than  the  DEM  data,  for  each  calculation,  DEM  height  values  are  found  by 

performing  a  two  dimensional  interpolation  of  the  DEM  data  at  the  locations  of  the 
shifted  and  rotated  SAR  pixel  locations  [143].  A  finer  spatial  registration  is  done  by 
shifting  the  position  of  {is,js)  by  offsets  smaller  that  the  DEM  pixel  resolutions  until 
a  minimum  is  found. 

To  adjust  for  orbit  errors  that  manifest  themselves  in  uniform  resolution  errors 
[155],  the  SAR  resolutions  in  both  the  range  and  azimuthal  directions  are  slightly 
modified  until  a  minimum  is  found.  The  resolutions  are  modified  with 


^XsAR  =  Ax, 


1 


SAR 


1  +  kx 


(5.57) 


and 


^USAR  =  ^VSAR 


1  -{•  ky 


(5.58) 


by  shifting  kx  and  ky  small  amounts  about  zero.  It  is  important  to  note  each  time  the 
range  resolution  is  modified  the  height  must  be  recalculated  since  the  height  inversion 
algorithm  depends  on  the  range  resolution. 

The  last  step  is  to  find  a  height  offset  that  can  be  subtracted  from  the  SAR 


172 


CHAPTER  5.  PHASE  UNWRAPPING  OF  SAR  INTERFEROMETRY 


1.  Invert  height  and  forshorten  data. 

2.  Provide  initial  SAR  range  and  azimuthal  resolutions. 

3.  Provide  initial  DEM/SAR  pixel  match. 

4.  Provide  initial  offset  angle. 

5.  Find  optimum  pixel  match  by  minimizing  RMS  height  error. 

6.  Find  optimum  angle  by  minimizing  RMS  height  error. 

7.  Find  optimum  sub-pixel  match  by  minimizing  RMS  height  error. 

8.  Find  optimum  resolution  scales. 

9.  Adjust  SAR  data  to  match  DEM  mean  height. 

Table  5.9:  Registration  steps. 


Figure  5-12:  Forshortening  diagram. 

inverted  height  so  that  the  mean  heights  of  both  data  subsets  are  equal.  A  summary 
of  registration  steps  is  listed  in  Table  5.9. 


5.9.3  Forshortening 

As  mentioned,  forshortening  adjusts  the  position  of  each  pixel  along  the  range 
direction  to  compensate  for  the  fact  that  range  pixels  are  not  necessarily  in  order  as 
far  as  the  ground  dimension  is  concerned.  Figure  5-12  is  a  diagram  that  contains 
the  necessary  elements  to  explain  the  forshortening  process.  The  horizontal  axis 
represents  the  range  direction  and  the  vertical  axis  represents  terrain  height.  The 
diagonal  lines  represent  four  range  bins,  labeled  1,  2,  3,  and  4,  at  a  constant  azimuthal 
step.  Any  scatterer  inside  a  range  bin  will  return  electromagnetic  energy  used  to 
calculate  the  height  for  that  specific  azimuthal  step  and  range  distance. 
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Consider  the  three  scatters,  A-C,  A  and  B  are  at  the  reference  height  of  zero  and 
C  is  at  a  height,  h  above  the  reference  height.  If  the  terrain  is  flat  and  includes 
scatters  A  and  B,  their  returns  will  be  in  range  bins  2  and  3  and  thus  converted  into 
the  proper  spatial  arrangement  i.e.  A  is  to  the  left  of  B.  However,  if  the  terrain  is  a 
hill  and  slanted  upward  including  scatters  A  and  C,  the  returns  will  be  in  range  bins 
2  and  1.  In  other  words,  scatterer  C  will  be  placed  to  the  left  of  A  because  it  was  in 
a  lower  numbered  range  bin.  The  result  is  a  mismatch  between  the  inverted  image 
map  and  the  true  map. 

The  spatial  errors  can  be  fixed  by  calculating  the  projection  of  the  inverted  height, 
g,  along  the  ground  using  * 

g  —  h  tan  •  (5.59) 

The  position  of  the  pixel  is  then  shifted  by  g.  If  g  is  positive,  the  pixel  height  is  shifted 
to  the  right  and  if  g  is  negative,  the  pixel  height  is  shifted  to  the  left.  In  this  way, 
all  pixels  are  shifted  in  the  image.  Equal  pixel  spacing  is  returned  by  interpolation 
at  the  ground  range  resolution  of  the  SAR  configuration. 


5.9.4  SAR/DEM  Registration 

Now  that  the  registration  process  has  been  described,  this  section  will  provide  the 
details  of  the  subset  of  data  used  to  compare  the  the  phase  unwrapping  techniques. 
The  registration  was  performed  using  a  300  x  200  set  of  DEM  data  points  taken  from 
the  northern  portion  of  DEM  image.  Figure  5-13  shows  the  relative  position  and 
Figure  5-14  is  a  gray  scale  representation  of  of  the  DEM  subset  with  white  and  black 
representing  the  highest  and  lowest  elevations  respectively.  The  SAR  subset  used  for 
registration  included  150  x  150  pixels  as  shown  in  Figure  5-15.  The  terrain  features 
can  easily  be  seen  in  both  images.  This  is  the  first  set  of  ERS  data  to  be  unwrapped; 
we  will  refer  to  it  as  ERS-1  Data  Set  1. 
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(1201) 


Figure  5-13:  DEM  subset:  shaded  square  indicates  subset  used  for  registration. 


(300) 


Figure  5-14:  Gray  scale  image  of  Alaskan  DEM  terrain  used  for  registration. 
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Figure  5-16:  Weighting  mask  based  on  coherence  map,  (ERS-1  Data  Set  1). 


5.9.5  Determination  of  data  weights 


The  three  types  of  weights  described  in  Sections  5.3.2  and  5.6  were  used  to  unwrap 
the  real  phase  data  with  the  weighted  least  squares  method.  Figures  5-16  through 
5-20  represent  the  weights  applied  to  the  interferogram  of  Figure  5-15.  In  these  gray 
scale  pictures  or  masks  white  represents  the  largest  weight  of  1,  black  represents  the 
smallest  weight  of  0,  and  the  grays  represent  continuous  weights  between  1  and  0. 

Figure  5-16  is  the  mask  based  on  straight  coherence  data  and  Figures  5-17,  5-18, 
and  5-19  are  the  modified  coherence  masks  with  coherence  thresholds,  Cthr,  of  0.8, 
0.5,  and  0.2.  Setting  Cthr  =  0  is  equivalent  to  no  weighting  {i.e.  a  straight  least 
square  method)  and  setting  Cthr  =  1  is  equivalent  to  straight  coherence  weighting 
(Figure  5-16). 
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Figure  5-17:  Weighting  mask  based  on  coherence  map  with  Cthr  —  0.8,  (ERS-1 
Data  Set  1). 


Figure  5-18:  Weighting  mask  based  on  coherence  map  with  Cthr  =  0.5,  (ERS-1 
Data  Set  1). 
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Figure  5-19:  Weighting  mask  based  on  coherence  map  with  Cthr  =  0.2.  (ERS-1 
Data  Set  1). 


Figure  5-20:  Weighting  mask  based  on  residues  (ERS-1  Data  Set  1). 
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5.9.6  Results 

Data  Set  1 

This  section  includes  a  compilation  of  the  height  errors  associated  with  four  specific 
phase  unwrapping  processes.  The  first  is  based  on  the  straight  least  squares  method, 
Figure  5-21;  the  second  and  third  are  based  on  the  weighted  least  squares  method 
using  the  coherence  weighting  and  residue  weighting  schemes,  Figures  5-23  and  5-27 
respectively;  and  the  last  is  based  on  the  global  optimal  branch  cut  method.  Figure  5- 
24.  The  gray  scale  plots  represent  the  errors  based  on  location  with  white  being  the 
largest  positive  error  of  125  m  and  black  being  the  largest  negative  error  of  -125  m. 

Associated  with  each  gray  scale  plot  is  a  histogram  of  the  pixel  height  errors,  Fig- 
mes  5-25,  5-26,  5-27,  and  5-28.  The  spikes  at  zero  are  an  artifact  of  the  last  step  of  the 
registration  process  (Table  5.9)  that  matches  the  mean  heights  of  the  SAR  and  DEM 
terrain  images.  On  ERS-1  Data  Set  1,  the  least  squares  methods  impose  a  global 
error  in  the  inverted  terrain  height  as  indicated  by  the  peak  in  the  error  to  the  left  of 
the  0  m  error  mark  on  the  horizontal  axis.  The  unweighted  least  squares  unwrapping 
process  produced  at  peak  at  40  m,  while  the  coherence  and  residue  weighting  schemes 
produced  peaks  at  25  m  and  8  m  respectively.  The  optimal  branch  cut  method  pro¬ 
duces  local  errors  indicated  by  the  symmetric  error  histogram  about  the  Om  error 
mark. 

Using  a  weighting  least  squares  method  reduces  the  global  error  of  the  unweighted 
least  squares  method.  The  new  residue  weighting  scheme  offered  much  better  improve¬ 
ment  over  the  weighting  based  on  coherence  data.  Figure  5-29  represents  the 

hrms  error  for  the  inverted  height  based  on  the  weighted  least  squares  method  using 
the  modified  coherence  weighting  scheme  on  ERS-1  Data  Set  1.  Over  a  twenty  per¬ 
cent  reduction  in  the  firms  error  was  realized  over  the  no- weight  least  squares  method 
(^Cthr  =  0  when  using  this  method.  Figure  5-30  compares  the  branch  cut  method  and 
the  weighted  least  squares  method  versus  iterations  of  the  Picard  iteration  method. 
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Figure  5-22: 
Data  Set  1). 


Least  Square  Plot 


Figure  5-21:  Least  square  error  plot  (ERS-1  Data  Set  1). 


Weighted  Least  Square  (Coherence)  Plot 


Weighted  least  square  error  plot  using  coherence  mask  (ERS-1 
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Weighted  Least  Square  (Residue)  Error  Plot 


Figure  5-23:  Weighted  least  square  error  plot  with  residue  mask  (ERS-1  Data  Set  1). 


Branch  Cut  Error  Plot 


Figure  5-24:  Branch  cut  error  plot  (ERS-1  Data  Set  1) 
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Figure  5-25:  Histogram  of  least  square  height  errors  (ERS-1  Data  Set  1).  Standard 
deviation:  38.60  m. 


Figure  5-26:  Histogram  of  weighted  least  square  height  errors  using  coherence  mask 
(ERS-1  Data  Set  1).  Standard  deviation:  29.47m. 
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Figure  5-27:  Histogram  of  weighted  least  square  height  errors  using  residue  mask 
(ERS-1  Data  Set  1).  Standard  deviation:  18.26m. 


Figure  5-28:  Histogram  of  branch  cut  height  errors  (ERS-1  Data  Set  1).  Standard 
deviation:  15.15  m. 


184 


CHAPTER  5.  PHASE  UNWRAPPING  OF  SAR  INTERFEROMETRY 


Figure  5-29:  Comparison  of  hrms  error  versus  Cthr  setting  of  height  inverted  with 
weighted  least  squares  with  coherence  weighting  on  ERS-1  Data  Set  1. 

On  Data  Set  1  data,  the  branch  cut  method  performed  exceptionally  well  with  an 
error  of  15.15  m.  The  straight  least  squares  method  has  an  error  of  38.60  m.  The 
weighted  least  squares  method  with  residue  weighting  has  an  error  of  18.26  m.  And 
finally,  the  weighted  least  squares  with  coherence  weighting  method  has  an  error  of 
29.47  m.  It  appears  that  the  branch  cut  method  and  the  hybrid  weighted  least  squares 
method  offered  the  best  performance. 

Data  Set  2 

In  this  section  another  set  of  ERS-1  data  is  use  to  study  the  various  unwrapping 
techniques.  We’ll  call  this  data  ERS-1  Data  Set  2.  The  interferogram  and  its  corre¬ 
sponding  coherence  data  are  shown  in  Figures  5-31  and  5-32.  This  data  set  contains 
256  X  256  pixels  and  the  entire  image  is  unwrapped.  The  corresponding  DEM  ter¬ 
rain  feature  can  be  seen  as  the  mountain  peaks  at  the  lower  center  of  the  map  in 
Figure  5-14.  The  residue  mask  is  represented  in  Figure  5-33.  Instead  of  showing 
the  grayscale  error  plots,  just  the  height  error  histograms  are  shown  here.  The  first 
histogram  represents  the  height  error  when  the  straight  least  squares  method  is  used. 


Figure  5-31:  ERS-1  Data  Set  2. 
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Figure  5-32:  Coherence  data  for  ERS-1  Data  Set  2. 


Figure  5-33;  Weighting  mask  based  on  residues  for  ERS-1  Data  Set  2. 
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Method 

(“) 

Time  (s) 

LS 

74.85 

8.43 

W-LS  (coherence) 

39.18 

43.64 

W-LS  (residues) 

20.91 

38.87 

Global  optimum  BC 

22.30 

588.98 

Local  optimum  BC 

22.40 

20.55 

Table  5.10:  Comparison  of  rms  height  error  and  CPU  times  for  the  phase  unwrapping 
methods.  LS:  least  squares,  W-LS:  weighted  least  squares  and  BC:  branch  cut. 


Figure  5-35;  the  second  and  third  histograms  are  based  on  the  weighted  least  squares 
method  using  the  coherence  weighting  and  residue  weighting  schemes,  Figures  5-36 
and  5-37  respectively.  The  fourth  histogram  is  based  on  the  global  optimal  branch  cut 
method,  Figure  5-38  and  the  last  histogram,  Figure  5-39,  represents  the  height  errors 
due  to  a  local  optimum  branch  cut  method.  Table  5.10  lists  the  standard  deviation, 
(KTrlr)^  for  each  method  used.  In  addition  to  the  error,  the  CPU  times  are  provided. 
These  times  are  meant  to  show  relative  times.  None  of  the  code  was  optimized  to 
provide  the  fastest  execution  times. 

The  results  of  unwrapping  ERS-1  Data  Set  2  show  that  the  least  squares  method 
is  not  always  appropriate  for  unwrapping  SAR  interferograms  the  rms  height  error 
is  almost  twice  the  next  largest  error.  The  weighted  least  squares  method  using 
the  coherence  data  weighting  scheme  offers  great  improvement  over  straight  least 
squares,  but  as  before,  the  residue  scheme  offers  an  still  an  additional  improvement. 
The  residue  weighted  least  squares  and  the  global  and  local  branch  cut  methods  offer 
the  best  rms  errors  of  20.91m,  22.30  m,  and  22.40  m  respectively.  However,  their 
respective  CPU  times  are  38.87  s,  588.98  s  and  20.55  s.  The  local  branch  cut  method 
is  the  fastest  followed  by  the  weighted  least  squares,  but  the  global  branch  cut  method 
is  very  slow.  The  data  set  contains  2663  residue  pairs  and  finding  the  optimum  takes 
almost  10  minutes.  Figure  5-34  are  grayscale  pictures  of  the  unwrapped  phase  of 
Data  Set  2  using  the  branch  cut  method  with  residue  connections  made  with  nearest 
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neighbor,  local  optimum,  and  global  optimum  branch  cut  algorithms.  One  can  clearly 
see  the  effects  of  long  branch  cuts  in  the  top  image  (nearest  neighbor).  Those  phase 
errors  are  removed  with  a  local  optimum  scheme  in  the  center  image  which  is  nearly 
identical  to  the  lower  image  (global).  The  center  image  was  unwrapped  in  less  than 
4  percent  the  time  necessary  to  unwrap  the  lower  image  with  the  global  algorithm. 

As  before  with  ERS-1  Data  Set  1,  using  a  weighting  least  squares  method  reduces 
the  global  error  of  the  unweighted  least  squares  method.  The  new  residue  weighting 
scheme  offered  considerable  improvement  over  the  weighting  based  on  coherence  data. 


Data  Set  3 

In  this  section  another  set  of  ERS-1  data  is  used  to  study  the  various  unwrapping 
techniques.  We’ll  call  this  data  ERS-1  Data  Set  3.  The  interferogram  and  its  corre¬ 
sponding  coherence  data  are  shown  in  Figures  5-40  and  5-41.  This  data  set  contains 
256  X  256  pixels  and  the  entire  image  is  unwrapped.  The  corresponding  DEM  terrain 
feature  can  be  seen  as  the  horse  shoe  shaped  mountain  peak  at  the  center  of  the 
map  in  Figure  5-14.  The  residue  mask  is  represented  in  Figure  5-42.  The  data 
was  selected  to  show  the  difficulty  even  the  best  methods  have  with  high  fringe  den¬ 
sity  compounded  with  low  coherence.  Data  Set  3  was  unwrapped  with  the  hybrid 
weighted  least  squares  method  using  the  residue  mask  in  Figure  5-42.  The  results  are 
shown  in  Figures  5-43  and  5-44.  The  unwrapping  lost  most  of  the  height  associated 
the  many  fringes  that  wrap  around  the  lower  portion  of  the  mountain  peak.  The 
result  is  that  the  peak  height  is  underestimated  by  over  120  m  and  the  flat  portion 
at  the  base  of  the  mountain  is  overestimated.  The  histogram  shows  a  global  error  at 
15  m.  The  global  optimum  branch  cut  method  offers  similar  rms  height  error,  39.07  m 
versus  39.33  m  for  the  hybrid  method.  Several  large  portions  of  the  map  are  over 
estimated  represented  by  the  light  patches  on  lower  portion  of  Figure  5-45.  Even 
with  the  addition  of  many  image  residues  on  the  boundary  were  not  able  to  overcome 


5.9.  APPLICATION  AND  ANALYSIS  OF  TECHNIQUES  ON  REAL  DATA  189 


Figure  5-34:  Comparison  of  nearest  neighbor  (top),  local  optimum  (center)  and  global 
optimum  (bottom)  branch  cut  phase  unwrapping  for  ERS-1  Data  Set  2. 
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Error  (m) 


Figure  5-35:  Height  error  histogram  for  straight  least  squares  method  applied  to 
ERS-1  Data  Set  2.  Standard  deviation:  74.85  m. 


Figure  5-36:  Height  error  histogram  for  weighted  least  squares  method  with  coherence 
weighting  applied  to  ERS-1  Data  Set  2.  Standard  deviation:  39.18  m. 
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Figure  5-37:  Height  error  histogram  for  weighted  least  squares  method  with  residue 
weighting  applied  to  ERS-1  Data  Set  2.  Standard  deviation:  20.91  m. 


Figure  5-38:  Height  error  histogram  for  global  optimum  branch  cut  method  applied 
to  ERS-1  Data  Set  2.  Standard  deviation:  22.30  m. 
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Figure  5-39:  Height  error  histogram  for  local  optimum  branch  cut  method  applied  to 
ERS-1  Data  Set  2.  Standard  deviation:  22.40  m. 


Figure  5-40:  ERS-1  Data  Set  3. 
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Figure  5-43:  Height  error  for  hybrid  weighted  least  squares  method  applied  to  ERS-1 
Data  Set  3.  Standard  deviation:  39.33  m. 

the  errors.  An  image  residue  was  placed  at  the  boundary  for  each  true  residue  within 
5  pixels  from  the  interferogram’s  edge. 


5.10  Summary 

In  this  chapter,  different  phase  unwrapping  schemes  were  described  and  applied  to 
both  simulated  and  real  SAR  interferograms.  The  techniques  used  were  weighted  and 
unweighted  least  squares  and  optimal  branch  cut  global  and  local  methods.  When 
no  noise  is  present  phase  unwrapping  is  straight  forward  and  uninteresting;  however, 
when  noise  is  present  the  phase  unwrapping  process  in  non-trivial  and  there  is  much 
interest  in  proper  treatment  of  the  noise. 

Although  the  least  squares  method  is  very  fast  because  of  the  Discrete  Cosine 
Transform  formulation,  it  does  not  treat  noisy  data  very  well  as  shown  when  at¬ 
tempting  to  unwrap  the  simulated  interferogram.  Once  the  rms  phase  noise  reached 
35°,  the  least  squares  method  was  unable  to  unwrap  the  phase  data  effectively.  Up 
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Figure  5-44;  Height  error  histogram  for  hybrid  weighted  least  squares  method  applied 
to  ERS-1  Data  Set  3.  Standard  deviation:  39.33  m. 


Figure  5-45:  Height  error  for  global  optimum  branch  cut  method  applied  to  ERS-1 
Data  Set  3.  Standard  deviation:  39.07  m. 
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Figure  5-46:  Height  error  histogram  for  global  optimum  branch  cut  method  applied 
to  ERS-1  Data  Set  3.  Standard  deviation:  39.07  m. 

to  35°  rms  phase  noise,  all  phase  unwrapping  algorithms  performed  equally  well  be¬ 
cause  the  noise  was  not  great  enough  to  produce  phase  inconsistencies  in  the  form 
of  residues.  At  35°  rms  noise  there  were  only  34  residues  in  the  image  of  65536  pix¬ 
els.  However  at  55°  rms  noise  there  were  4983  residues  in  the  image  of  65536  pixels 
and  both  the  residue  (hybrid)  weighted  least  squares  and  the  optimum  branch  cut 
methods  kept  the  rms  phase  error  to  less  than  one  radian. 

Both  the  hybrid  weighted  least  squares  and  the  optimum  branch  cut  methods  of¬ 
fered  essentially  the  same  error  performance  when  applied  to  both  simulated  SAR  in- 
terferograms  with  additive  noise  and  real  ERS-1  interferograms.  These  methods  con¬ 
sistently  showed  improvement  over  the  straight  least  squares  and  coherence  weighted 
least  squares  methods.  Specifically,  the  methods  offer  over  fifty  percent  reduction  in 
root  mean  square  (rms)  height  error  compared  to  the  straight  least  squares  method 
and  over  thirty-five  percent  reduction  in  rms  height  error  compared  to  the  weighted 
least  squares  method  based  on  coherence  data  weighting  schemes. 

The  global  optimum  branch  cut  method  is  appropriate  when  there  are  not  that 
many  residues;  however,  the  local  optimum  provide  near  identical  results  in  far  less 
time.  The  branch  cut  methods  is  more  likely  to  produce  local  errors  as  shown  by  the 
symmetrical  error  histograms  and  the  hybrid  weighted  least  squares  method  produces 
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a  fairly  symmetrical  histogram  especially  when  compared  to  the  straight  least  squares 
and  coherence  weighted  least  squares  methods.  In  other  words,  using  a  weighting 
scheme  greatly  reduces  the  global  error  introduced  by  the  least  squares  method. 

The  hybrid  method  could  be  improved  by  using  a  different  algorithm  to  solve  the 
weighted  least  squares  problem  that  guarantees  convergence.  Although  the  Picard 
iteration  always  converged  in  this  application,  simply  selecting  an  iteration  number 
is  not  the  best  way  to  solve  the  problem.  A  conjugate-gradient  solver  as  proposed  in 
[170]  would  guarantee  convergence  and  prevent  running  unnecessary  iterations. 

The  branch  cut  method  could  be  improved  by  characterizing  residues  by  their 
source.  For  example,  residue  pairs  due  to  speckle  noise  are  close  to  each  other  and  han¬ 
dled  quite  effectively  with  the  optimum  methods.  However,  residues  created  by  under 
sampled  fringes  caused  by  very  steep  terrain  features  can  produce  residue  pairs  rela¬ 
tively  far  apart  and  are  not  treated  properly  in  the  presence  of  many  other  residues. 
If  these  residues  can  be  isolated  they  can  be  handled  separately  as  proposed  in  [163]. 
One  way  these  residues  could  be  found  is  by  using  the  local  optimum  branch  cut 
method.  When  the  local  optimum  algorithm  cannot  reduce  the  branch  cut  length  of 
a  long  cut,  it  probably  means  there  is  a  non-noise  produced  residue  in  the  vicinity  of 
the  long  cut.  These  areas  can  be  flagged  and  treated  with  other  techniques. 
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Chapter  6 
Conclusion 


In  this  thesis,  applications  of  numerical  techniques  to  electromagnetic  problems  in 
microelectronic  and  radar  imaging  systems  were  investigated.  Dielectric  rib  waveg¬ 
uide  discontinuities  were  analyzed  with  the  Finite  Difference  Time  Domain  (FDTD) 
method.  Within  the  FDTD  framework,  two  techniques  were  applied  to  the  study  of 
dielectric  rib  waveguides.  First,  Berenger’s  Perfectly  Matched  Layer  was  used  to  trun¬ 
cate  the  multi-layered  dielectric  structure.  Second,  the  rib  waveguide’s  fundamental 
mode  was  determined  with  a  two  dimensional  FDTD  simulation  for  subsequent  use 
in  a  three  dimensional  simulation.  These  numerical  techniques  were  used  to  study 
the  effect  of  bend  discontinuities  in  the  rib  waveguide. 

An  Impedance  Boundary  Condition  (IBC)  was  developed  for  two  dimensional 
FDTD  simulations  that  can  replace  thin  sheets  of  highly  conducting  material.  The 
IBC  accurately  models  the  conductor  loss  over  a  wide  frequency  range.  In  this  way,  the 
FDTD  method  is  improved  by  allowing  a  method  to  model  conductor  loss  without  the 
need  for  very  large  computational  domains.  The  two  dimensional  IBC  was  validated 
through  the  comparison  of  resonant  cavity  quality  factors  determined  with  the  FDTD 
simulation  data  and  those  calculated  with  analytical  methods. 

Phase  unwrapping  techniques  for  the  inversion  of  terrain  height  using  Synthetic 
Aperture  Radar  Interferometry  (InSAR)  data  were  analyzed.  A  hybrid  phase  unwrap- 
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ping  technique  that  combines  the  weighted  least  squares  method  with  the  residues  of 
the  branch  cut  method  was  presented.  Optimal  branch  cut  methods  were  also  pre¬ 
sented:  a  global  one  that  includes  all  residues  before  placing  any  branch  cuts  and  a 
local  one  that  uses  a  subset  of  residues  to  fix  long  branch  cuts  introduced  when  using 
a  nearest  neighbor  connection  scheme.  These  new  SAR  phase  unwrapping  methods 
were  used  to  unwrap  simulated  and  real  SAR  interferograms  and  then  were  compared 
to  straight  least  squares  and  other  weighted  least  squares  methods. 


6.1  Analysis  of  Dielectric  Waveguide  Discontinu¬ 
ities 

A  study  of  the  Perfectly  Matched  Layer  showed  that  the  PML  conductivities  used  to 
absorb  the  outgoing  waves  in  a  multiple  dielectric  layer  simulation  must  not  be  as¬ 
signed  independently  with  the  method  of  normal  incidence  used  by  Berenger.  Instead, 
once  one  of  the  layer’s  conductivities  is  calculated  the  others  must  be  constructed  so 
that  the  loss  tangents  in  each  layer  are  equal.  With  this  done,  there  will  be  no  re¬ 
flection  coefficient  singularities  in  both  the  normal  and  tangential  directions  and  the 
PML  will  simulate  open  space  to  the  maximum  extent  possible  under  the  discretiza¬ 
tion  scheme.  If  the  interdependency  of  the  dielectric  layers  is  ignored  even  though 
each  layer  is  perfectly  matched,  the  material  differences  at  the  interfaces  act  as  sources 
for  unwanted  solutions  that  contaminate  the  simulation  volume.  With  the  matching 
conditions  determined  for  multi-layered  dielectrics,  dielectric  rib  waveguide  structmes 
were  terminated  with  the  Berenger  Perfectly  Matched  Layer  and  each  dielectric  layer 
was  matched  considering  all  the  layers  within  the  simulation  volume. 

It  was  shown  that  the  computational  domain  can  be  reduced  when  simulating 
dielectric  rib  waveguide  with  the  FDTD  numerical  technique  when  the  waveguide 
fundamental  mode’s  spatial  distribution  is  calculated  first  and  used  to  excite  the 
guide.  The  better  excitation  means  that  shorter  distances  between  the  excitation 
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plane  and  the  discontinuity  are  needed  for  the  mode  to  settle.  A  reduction  of  twenty 
percent  was  shown.  This  reduction  implies  bigger  problems  can  be  simulated  within 
the  same  time  as  those  simulations  using  other  excitation  techniques,  or  the  same 
simulations  can  be  run  in  shorter  times. 

Although  the  dielectric  rib  waveguide  provides  very  low  loss  for  high  frequency 
signals  compared  to  microstrip  or  coplanar  structures,  the  rib  waveguide  bend  does 
introduce  a  significant  amount  of  loss.  As  the  abruptness  of  the  turn  is  lessened,  the 
rib  guides  more  energy  around  turn;  however,  gradual  turns  require  more  space  on 
an  integrated  circuit.  With  the  improvements  to  the  FDTD  method  presented,  rib 
bend  structures  can  be  studied  to  improve  their  loss  characteristics. 


6.2  Modeling  of  Thin  Finite  Conductivity  Sheets 

A  new  Impedance  Boundary  Condition  for  good  conductors,  satisfying  1  <<  was 
extended  for  use  in  two  dimensional  FDTD  simulations.  The  IBC  was  described  in 
great  detail  and  studied  in  one  dimensional  simulations.  The  influence  of  FDTD 
parameters  and  the  number  of  IBC  expansion  terms  on  the  accuracy  of  the  IBC  was 
presented.  The  results  showed  that  the  time  step  constraint  is  more  stringent  than 
the  Courant  stability  limit.  The  time  step  must  satisfy  At  < 

The  conductivity  range  studied  spanned  four  orders  of  magnitude  from  10“*  to 
Two  variations  of  the  IBC  were  investigated;  one  that  approximates 
the  convolution  with  a  piecewise  constant  magnetic  field  and  another  that  uses  a 
piecewise  linear  magnetic  field.  When  a  piecewise  constant  approximation  is  used 
for  the  magnetic  field,  the  amount  of  computer  resources  (as  measured  by  number  of 
expansion  terms)  is  directly  proportional  to  the  conductivity.  Lower  conductivities 
result  in  larger  pole  magnitudes  and  lesser  sensitivity  to  an  increased  number  of  terms. 
The  piecewise  constant  approximation  produces  an  error  that  cannot  be  overcome 
by  more  expansion  terms.  The  errors  from  the  impedance  approximation  and  the 
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magnetic  field  approximation  combine  to  form  an  optimum  number  of  expansion 
terms.  Using  more  than  the  optimum  number  of  terms  will  result  in  less  accuracy. 
On  the  other  hand,  with  a  piecewise  linear  approximation,  the  error  to  the  magnetic 
field  is  greatly  reduced  and  an  optimum  number  of  terms  is  not  evident.  The  more 
terms  used  the  more  accurate  the  IBC;  however,  the  error  quickly  levels  off. 

When  extended  to  two  dimensions,  the  impedance  boundary  condition  works  very 
well  over  the  wide  range  of  skin  depths  studied,  3  185.  The  lowest  conductivity, 

a  =  5.8  X  10^  offered  the  biggest  error  of  almost  13%,  but  the  error  was  very 

small  at  conductivities  near  those  of  the  most  popular  conductors  of  silver,  copper, 
gold  and  aluminum.  The  IBC  is  an  excellent  model  for  thin  sheets  of  good  conductors 
with  finite  conductivity.  The  IBC  offers  the  advantage  of  working  over  a  wide  range 
of  frequencies  with  very  little  computational  overhead.  For  typical  conductors,  like 
copper,  125  terms  per  tangential  electric  field  component  are  needed  with  errors  in 
Q  less  than  5%. 


6.3  Phase  Unwrapping  of  SAR  Interferograms 

Numerical  techniques  to  unwrap  the  phase  of  synthetic  aperture  radar  interferograms 
were  investigated.  The  foundations  of  phase  unwrapping  were  described,  specifically 
least  squares  and  optimal  branch  cut  unwrapping.  New  modifications  to  these  tech¬ 
niques  were  applied  to  both  simulated  and  real  interferograms.  After  the  application 
of  these  techniques  to  the  real  SAR  data,  the  phase  was  used  to  invert  terrain  height 
and  the  height  was  compared  to  ground  truth. 

Although  the  least  squares  method  is  very  fast  because  of  the  Discrete  Cosine 
Transform  formulation,  it  does  not  treat  noisy  data  very  well.  In  simulated  SAR 
interferograms,  once  the  rms  phase  noise  reached  35°,  the  least  squares  method  was 
unable  to  unwrap  the  phase  data  effectively.  Both  the  hybrid  weighted  least  squares 
and  the  optimum  braneh  cut  methods  offered  essentially  the  same  error  performance 
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when  applied  to  both  simulated  SAR  interferograms  with  additive  noise  and  real 
ERS-1  interferograms.  These  methods  consistently  showed  improvement  over  the 
straight  least  squares  and  coherence  weighted  least  squares  methods.  Specifically, 
the  methods  offer  over  fifty  percent  reduction  in  root  mean  square  (rms)  height  error 
compared  to  the  straight  least  squares  method  and  over  thirty-five  percent  reduction 
in  rms  height  error  compared  to  the  weighted  least  squares  method  based  on  coherence 
data  weighting  schemes. 

The  global  optimum  branch  cut  method  is  appropriate  when  there  are  not  that 
many  residues;  however,  the  local  optimum  method  provides  near  identical  results 
in  far  less  time.  The  branch  cut  methods  are  more  likely  to  produce  local  errors  as 
shown  by  the  symmetrical  error  histograms  and  the  hybrid  weighted  least  squares 
method  produces  a  fairly  symmetrical  histogram  especially  when  compared  to  the 
straight  least  squares  and  coherence  weighted  least  squares  methods.  In  other  words, 
using  a  weighting  scheme  greatly  reduces  the  global  error  introduced  by  the  least 
squares  method.  In  fact,  the  hybrid  method  closely  approaches  the  desirable  local 
error  distribution  of  the  branch  cut  method. 


6.4  Outlook 

Although  much  was  accomplished  in  this  work,  improvements  can  always  be  made 
and  more  could  be  done.  The  most  logical  extension  of  the  dielectric  waveguides 
work  would  be  to  thoroughly  investigate  ways  to  reduce  the  loss  through  the  turn 
discontinuity.  Different  geometries  and  materials  could  be  tested  to  reduce  the  loss 
while  keeping  the  size  small.  In  the  Impedance  Boundary  Condition  work,  the  IBC 
could  be  extended  to  three  dimensions.  Finally,  the  InSAR  weighted  least  squares 
phase  unwrapping  technique,  in  this  work,  could  be  improved  by  replacing  the  Picard 
iteration  algorithm  with  another  solver  that  will  improve  the  convergence  properties. 
The  optimal  branch  cut  methods  could  be  improved  by  characterizing  the  residues  by 
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there  most  likely  source.  In  this  way,  residue  pairs  that  may  be  relatively  far  apart 
can  be  handled  separately  leaving  the  optimum  algorithm  to  connect  the  pairs  that 
should  be  close  together. 

Despite  the  increases  in  computational  power  that  result  from  improving  computer 
technology,  electromagnetic  problems  require  large  amounts  of  computer  resources. 
Finding  better  numerical  techniques  or  improving  existing  numerical  techniques  will 
continue  to  be  worthy  goals. 


Appendix  A 


Mur  FDTD  Equations 


Given  a  computational  domain  with  boundaries  x  =  0,Nx,  y  =  0,Ny  and  z  =  0,Nz, 
the  tangential  electric  fields  at  all  six  boundary  faces  must  be  updated  with  the  ABC 
equations.  For  completeness,  all  the  field  equations  are  listed  below. 


A.l  1st  Order  Mur  FDTD  Equations 


+  l,N„k)  =  E;{i  +  l,N,-l,k) 

+  ^^lE;+'{i  +  \.N,-l,k)-E",{i  +  ^N„k)],  (A,2) 
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£^r'(i  +  i.i  N.)  =  BJCi  +  5, i  w,  -  1) 

'''  CoA(  + (>+ -  1) -BJ(x+ 


{A.4) 


(A.5) 


+  i  «:)  =  £;;(jv.  -  l,i  +  i  *) 


+  + 5’ +  (a.6) 


isr'('.:'  +  5’0)  =  -B,"(<.:'‘  +  5.i) 


(A.7) 


+  5.  JVJ  =  E^{i,j  +  5.  JV.  -  1) 


;.JV.)1, 


(A.8) 


i?^‘(0,i,*;  +  i)  =  £?(!, i,i:+i) 


(A.9) 


K*\N„j,  *  +  i)  =  Er(Af,  -  IJ,  A;  +  i) 

+  +  5)  -  +  5>1’ 


(A.10) 
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(A.11) 


£;,“+■(*  +  i,  Ny,k)  =  E^(i  +  \,Ny-  1,  k) 

+  +  i. TV,  -  u)  -  i?(i  +  i,  JV.,  m. 


(A.12) 


A. 2  2nd  Order  Mur  FDTD  Equations 

The  second  order  Mur  equations  are  substantially  more  complicated  than  the  first 
order  Mur,  they  are 


ET'(i  +  J.i.O)  =  +  5.  j.  1)  (A13) 

+£:(i + 1)  -  2£j(i + i,  j.  1) + £;(>  -  5,  j,  1)1 

+  2(A.tw  ^  ^  ^  °> 

+£;(•■  +  5,2  +  1, 1)  -  2£:(i  +  i  2, 1)  +  £?(•  +  i,  J  -  1. 1)1. 


£r'(i + 5.  j.  N.)  =  + 5,2,  Ny  - 1) 

9  A  r  1  1 

+  ^JA?T^'^-(*  2'^-  +  2'^'-  -  ')! 


(A.  14) 
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+E;(i  +  |i,  AT,  -  1)  -  2E:(i  +  i  J.  JV,  -  1)  +  E;(i  -  i  j,  N,  -  1)] 

+  2(Apn^t!+AU^(‘  +  5’^'  +  +  h’  - 

+£? (i  +  j. J  +  1.  JV.  -  1)  -  2£;(i  +  ij,  JV.  -  1)  +  £?(i  +  i. j  -  1,  Af,  -  1)1, 


£.”+>(»  +  J.  0.  k)  =  -£;-‘(i  +  5,1,4)  (A.15) 

CoAt  +  +  2’  *)  +  ^<'  +  2’ 

'*'  2(Aa:)2(coAl  + Aj/)^"®"^'  '*'  2’ ■•■  J’ “  2’ 

+£?(!  +  1 1, 4)  -  2B;(i  +  i,  1, 4)  +  £5(1  -  i,  1, 4)1 

'*'  2(A2)2(coAt  +  Asi)^^"*'^  2’*'’*"'''  2’°’*^'"''^^*^  2’*'’*^“ 

+£>”(>  +  5. 1, 4  +  1)  -  2E:(i  +  5, 1, 4)  +  E;(i  +  i,  1, 4  -  1)1, 


E:*Hi  +  5,  JV„,  4)  =  -E^-'ii  +  5,  iV,  -  1, 4) 


(A.16) 


cpAt  -  Ay 
cpAt  +  Ay 
2Ay  , 


+  +o,Ny-  1,  k)  +  E^ii  +  -,Ny-  1,  k)] 


+  1-^)1 

2{Axy{coAt  +  Ay)  ^  2’  ~  2’  “  2’ 

+E:{i  +  l,Ny-  1,  fc)  -  2E:{i  +  l,Ny-  1,  k)  +  E^ii  -\,Ny-  1,  k)] 

2{Azy{coAt  +  Ay)  2’  ^  +  2’  +  2’  ^  ~ 

+E-{i  +  1, 1,  A:  +  1)  -  2£;,"(i  +  i  1,  k)  +  E^^i  +  \,l,k-  1)], 
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j  +  -,k)  —  Ey  I  2>'v  v - / 

^  2(AwL)[^(°'  ;  +  i  +  1)  -  2^(0, 2  +  i  +  i  ^  -  1) 

+E;{1,  j  +  ^,k  +  l)-  2E'^{l,j  +  l,k)  +  E^iU  j  +  l,k~l)] 

^  wt°A/+ f 4*)  + 

+£;(1.  j  +  2  )d)  -  2B;(1,j  +  i  *:)  +  £,"(!,:;■- i  m, 


{A.17) 


£r‘(^-i  +  i  *=)  =  -•Ey’HAf.  -  1,7  +  5,  *;) 

+  ^  5’*>' 

+  +  5'  *=)  +  i 

+  +  i  i 

+£;;(iv,  - 1,  j  +  ^,  fc  + 1)  -  - 1,  j  i  A:)  +  £;;;(iVx  - 1, i  +  ^,  ^  - 1 

+  2(A,twAx)l^?<^--^'  +  +  5' '=)  +  i '^> 

+£;;(iv,  -  i,i  +  I  i)  -  2£^(A',  - 1.7  +  i  fc)  +  e;(n,  - 1,7  -  5.  m, 


(AIS) 


i,(;)  +  EJ(Af„i  +  i  fc-1) 
(JV.-1.7  +  i,i-l)l 


+B^(N.-l,j  +  l,k)-2E^{N, 


E;*'{i,j  +  5,0)  =  -E^-Hi.i  +  5, 1) 


{A.19) 


■1,7  +  5, 


210 


APPENDIX  A.  MUR  FDTD  EQUATIONS 


+£,"(!  +  1,3 +  1,1)-  +  1  1)  +  £"(,  -  1,  J  +  1  1)] 

+  +  I®)  -  +  5’®)  +  -  i®) 

+E;(i,j  + 1, 1)  -  2E;(i,j  +  i,  1)  +  i  1)1, 


>  - 1,1  + 


+E,"(i, 


j  +  1 1)  -  2E;(i,j  +  i,  1)  +  E;(i,j  -  i. 


N.)  =  -Br'fei  +  J.  JV.  -  1)  ( 

+  - 1) + i  -  D) 

2^^  1  1 

Az(coAt)^  1  ^  ,1 

+  2(Ax)1(coa4az)‘^»('  +  +  2'  +  2' +  ) 

+E;{i  +  1,3  +  i  iV.  -  I)  -  2E;(i,3  +  i  AT,  -  1)  +  e;{i  -1,3  +  i  iV.  -  1)J 

A^fcoAt)^  r  r,,  3  ,  1  1 

+  -  2- 

+£J(i.i  + 1 JV.  - 1)  -  + 1  iv,  - 1)  +  E;{i,3  - 1,  w,  - 1)], 


2^^  1  1 
+  +  2’  +  2’ 

+  2(amZ  az)  m  +  +  2E;(i,3 

PE^ii  +  l,j  +  i  iV,  -  1)  -  2E^iiJ  +  i  iV,  -  1)  + 

Az(coAtY  r  r,,  3  ,  1 

+  ^A,)^(coAt  +  A.)  +  2’ +  2 

- 1)  -  2£;”(i,i  + 1  iv^  - 1)  +  R-( 


(A.20) 


:i-l,j  +  i  TV.) 


£;^+40,j,A:  + ^)  = -£;^  ^{l,j,k  +  ^)  (A.21) 

coAt  +  ^  ^ 

2Ax  1  1 

+  2{A^!^T^K(«'^-  <=  +  5)  -  2^(®.1.*=  +  5)  +  ^(®'^-  *  -  5) 

+E:(1,3,  *  +  |)  -  2£J(1.;,  A:  +  i)  +  £?(1,  j,  A:  -  i)] 

AxfcnAi)^  11  1 

+  W(U  +  A.)t^?(°-^'  +  1,*  +  5)  -  2B?(0,,'.  A;  +  -)  +  £."(0,,-  -  1,  *  +  5) 

+E?{1,1  +  1,  A:  +  i)  -  2B;(1.  j,  A.  +  i)  +  B;(1,  j  -  1,  A:  +  i)]. 
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t  +  1)  =  -  1,  j,  *:  +  i)  (A.22) 

+  +  5)  +  +  5>’ 

+  *=  +  5)  +  ^  5>1 

+  2(A.)WA.)I^"(^-^''  ^  +  5)  -  +  5)  +  ^  -  5) 

+E^{N,  -  l,j,  k  +  ^)-  2E^{N,  -  l,j,  k+^)  +  E^{N^  -  IJ,  k  -  ^)] 

+  2(A,twA.)‘^"<^-^'  +  ^  5>  -  ^  5)  +  -  1' *=  +  5) 

+£;”(Ar,  -  IJ,  k  +  ^)-  2E:{N,  -  1, j,  k  +  ^)  +  E^{N,  -lj-l,k  +  i)], 


0,k  +  ^)  =  -Er\i,  1>  ^  +  5) 

9Ai/  1 


(A.23) 


^+5)) 


'  ““c7 

+£;,"(i,  1,  i  + 1)  -  2£;,"(i,  1,  fc  +  5)  +  E;{i,  1,  fc  -  i)l 
+  2(Ax^t  a/+  a„)  ^  ^  +  5)  -  0. + 1)  +  g;(^  - 1, 0. .  + 1 1 

+£?(»  +  1, 1,  *;  +  i)  -  2£”(».  1,  *:  +  i)  +  £?(!  -  1. 1.  *:  +  j)!. 


£?+■(»,  Af».  *  +  5)  =  -Br‘(i.  Af,  -  1,  *:  +  i) 

+  M  +  5)  +  (•.  -  U  +  i)l 

+  +  5)  +  +  5)' 


(A.24) 
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+£?(>,  W,  -  1,  +  ?)  -  2E?(i,  AT,  -  1,  i  +  i)  +  i5(i,  JV,  -  1,  i  -  1)] 

+  2(AmZav)^^^‘  ^  '■  '^■'’ <^  +  5)  -  +  5)  +  £?(>  -  1.  at..  A;  +  i) 

+£;r(i  +  1,  iV,  -  1,  A:  +  i)  -  2£;,"(^,  iV,  -  1,  A;  +  i)  +  -  1,  AT,  -  1,  A:  +  i)]. 


Appendix  B 

FDTD  PML  Difference  Equations 


The  the  difference  equations  for  the  twelve  PML  fields  are  listed  below. 


(B.l) 


+  Cfii  +  \,j,  k)  [HT'Hi  +  i  j  +  l,k)  -  Hr  Hi  +  \3-\k 


Eit\i  +  7:,hk)  = 


Cr{i  +  l,}.k)EUi  +  l,j,k) 


(B.2) 


Cr(i  +  \,3,k)  +  \,3,k  +  i)  -  Hr'‘(i  +  \,3,k-  5) 


Er\i,3  +  l^,k) 


Cr{i.3  +  \,k)E;ji,3  +  \,k) 


(B.3) 


cr(*>2  +  k)  Hr^(i  +  \,3  +  5.  k)  -  -^,3  +  ^,k 


ErKi,3  +  ^,k)  =  Cni.3  +  l,k)E^Ai,3  +  l<k) 


(B.4) 


+  C'J'(!,y+ i,(:)  Hr^iJ  +  ^,k+^)  -  Hr^'iiJ  +  ^,k  -  - 
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ET{i,j,h+h  =  CT{i,3,k  +  hEi,{i,j,k+h 


+  Cr{i,j,  k  +  i)  +  \.i,  k+\)-  -\,hk+  i)' 


E:;\i,i,k  +  -) 


Cr(i,hk  +  \)E%(i,i,kN-) 


(B.6) 


ai’(i,j,k  +  \)  j  +  5.*:  +  i)  -  -  i  ft  +  1)‘ 


where  the  two  constants,  Cf®  and  C2  ®  a  —  x,y,  z,  are  given  by 


(B.7) 


Cr{i,j,k)  = 


il-Cr{i,Tk)) 
a%{i,  j,k)  Aa 


Similarly,  the  the  difference  equation  forthe  magnetic  fields  used  in  the  PML  region 
of  the  computational  domain  are 


ET”{h3  +  2’k'k  2'iklxy  ^{i,j  +  2’ ^ 2^ 

Cr(i,3  +  \’k+\)  [^(i,  j  +  l.i  +  i)  -  £."(».  *:  +  i) 


2 +  -,*;+ -)  =  C^"(i,3 +  ^,k  +  ^)Hxz ‘(i,j +  ^,k  +  (B.IO) 

+  CriiJ  +  \,k  +  \)  [E^(i,3  +  i  *  +  1)  -  E;{i,j  +  i,t)' 


H^Hi+l,3,k  +  ^)  =  +  +  +  i  i,fc  +  i)  (BJl) 

+  Cr(i  +  5,  j,  A  +  5)  [iS?(i  +  1,2,  *:  +  i)  -  E-{i,j,  k  +  1)' 
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Cr(>  +  \,h  k  +  +  \,hk  +  \)  (B.12) 


Cr(i  +  i,j,*^  +  i) 


K(i  +  \^},k+l)-BZ(i  +  ]^,i,k) 


■n-zx 


r-  1  ■  1 

(*+  2^3  +  2’^^ 


Ci'^{i  +  -,j  +  -,  k)Hzx  ^{'i'+  2^j  +  2>  (B.13) 

+  -,i  +  -, k)  \Ey{i  +  l,j  +  -,k)  -  Ey{i,j  +  k) 


(B.14) 


+  Cf”(i  +  i.:i  +  i  i) 


E;(i  +  ^j  +  i,k)-E;(i  +  ^,},k) 


where  the  two  constants,  C'f’”  and  are  given  by 


(B.15) 


and 


Cr{iJ,k) 


{1-Crii,3,k)) 

(T^{iJ,k)Ao: 


(B.16) 


The  PML  coefficients  are  a  function  of  position,  direction  and  the  field  type.  They 
are  functions  of  position  because  both  the  permittivity  and  the  PML  conductivity 
are  functions  of  position.  They  are  functions  of  direction  since  the  spatial  deriva¬ 
tive  approximations  are  functions  of  grid  size  which  in  general  are  different  for  each 
direction.  Finally,  they  are  functions  of  field  type  since  electric  and  magnetic  field 
equations  depend  on  electric  and  magnetic  conductivites  respectively. 

As  previously  stated  the  PML  conductivities  are  derived  by  the  continuous  equa¬ 
tion. 


^e,m  _  e,m 
'-'a  '^amax 


a  =  x,y,  z. 


(B.17) 
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where  aamax  is  the  maximum  conductivity  and  d  is  the  PML  thickness  such  that 
d  =  niayers^oi  [uiayers  IS  the  number  of  discretized  PML  layers).  The  discretized 
electric  conductivity  (in  the  y  direction  for  example)  is  calculated  by 

(B.18) 

Since  the  magnetic  conductivities  are  shifted  by  half  a  grid  space  compared  to  the 
electric  fields,  the  magnetic  conductivites  are  also  shifted  so  that 

j  +  -j  Ay)  =  cr^max  (j) 

As  equations  (B.18)  and  (B.19)  show,  the  discretized  conductivities  are  calculated  as 
the  average  conductivity  over  one  grid  dimension  centered  at  the  field  location. 

With  a  uniform  grid,  the  coefficient  direction  dependency  is  removed.  In  a  single 
media  environment,  the  coefficents  are  just  a  function  of  the  PML  layer  and  the  field 
type.  For  example,  Cf®,  Cf”*,  and  Cf"*  become 


Cf{l)  = 


cm 


(1  -  cm) 

<^max9il)^  ’ 


C^{1)  = 


and 


cr(o  = 


(1  -  crm 

<a»S{0A  ’ 


where  I  is  the  index  of  the  PML  layer  and  g{l)  is  a  grid  factor  defined  as 


da,  a  =  x,y,  orz. 


(B.20) 

(B.21) 

(B.22) 

(B.23) 

(B.24) 
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Since  one  of  the  PML  matching  conditions  is 

a’”  O'® 

fJ^o  eo 

the  first  coefficient  is  field  independent,  i.e.  Cf  —  —  Ci. 


(B.25) 


Furthermore,  if  there  is  another  media  with  (ei,^o)  then 


Cf2(0  = 


(1  -  cum 


Since  the  PML  matching  conditions  between  two  dielectrics  is 


Cl  Co 

we  see  that  the  first  coefficient  is  media  independent,  i.e. 


(B.26) 

(B.27) 


(B.28) 


Now  let  there  be  N  different  dielectrics  defined  by  their  permittivities,  e-y,7  = 
0, 1, 2, . . . ,  N’  using  freespace  as  a  baseline  {a  =  0)  then  the  final  equations  for  the 
formulation  of  FDTD  PML  coefficients  are 


C.yl{l)  =  Coiil) 


(B.29) 


and 


C^l)  = 


(l-C'oi(O) 


The  magnetic  coefficients  are  defined  by 


(B.30) 


C'.IW  = 


(l-goi(O) 


(B.31) 
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which  simplifies  to 


C5(i)  = 


(l-CoiW) 


(B.32) 


becoming  media  independent.  Note  that  when  —  0  then  the  PML  explicit  expo¬ 
nential  equations  become  standard  Maxwell’s  difference  equations,  where 


and 


Coi  =  1, 

(B.33) 

re 

~  e,A’ 

(B.34) 

At 

(B.35) 

Appendix  C 

Transportation  Problem 

The  transportation  problem  is  used  to  construct  a  solution  to  the  optimal  set  of 
branch  cuts  that  minimizes  the  total  branch  cut  length. 

C.l  Introduction 

The  transport  problem  is  a  specific  example  of  a  wide  class  of  linear  programming 
problems  [181],  Essentially  the  problem  is  to  minimize  the  cost  of  transporting  goods 
to  m  locations  with  n  supply  vehicles.  Mathematically  the  transportation  problem 
minimizes  the  total  cost  defined  by 

m  n 

2  =  (C.l) 

i— 1  i— 1 

within  the  constraint  that  a  supplier  (supply  node),  i,  can  supply  no  more  than  a* 
units, 

n 

^  ^  ^i,j  —  ^  I5  2,  .  .  .  ,  TTlj  (C.2) 

J=1 

and  each  receiver  (demand  node)  j  must  receive  at  least  bj  units, 

m 

i  =  (0.3) 

i=l 
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Figure  C-1:  Transportation  network. 


and  no  units  are  shipped  from  demand  nodes  to  supply  nodes, 

Xij>0,  i  =  j  =  l,2,...,n  (C.4) 

and  the  cost  of  shipping  one  unit  from  the  frh  supply  node  to  the  jth  demand  node 
is  Cij.  The  values  ai,  bj,  and  Cij,  are  nonnegative  integers. 

This  problem  can  be  constructed  as  a  directed  network  with  a  set  of  nodes,  V, 
and  a  set  of  edges,  E.  A  directed  edge  is  a  connection  between  two  nodes  that  has 
a  defined  flow  direction.  The  network  is  bipartite  and  complete.  Bipartite  means 
that  the  nodes  can  be  divided  into  two  distinct  groups;  in  this  case  the  supply  nodes 
(i  =  1,2,...,  m)  and  the  demand  nodes  (j  =  1,2,...,  n).  Completeness  implies 
all  the  nodes  are  somehow  connected  to  each  other.  In  the  transport  problem  each 
supply  node  has  an  edge  to  each  of  the  n  demand  nodes  and  each  edge  has  a  unit 
cost,  Cij,  associated  with  using  it.  So  the  whole  problem  boils  down  to  finding  the 
shipping  amounts,  a;,^,  that  minimize  the  total  cost  2.  Figure  C-1  is  a  graph  of  the 
transportation  problem. 


C.2.  MAXIMUM  FLOW  THROUGH  A  NETWORK 
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The  problem  can  be  solved  only  if  there  are  enough  units  in  supply  to  meet  the 
demand, 

m  n 

i=i  j=i 

where  the  bj  is  the  demand  at  node  j.  All  the  deliveries  to  demand  node  j  should  add 
up  to  the  total  demand, 

n 

i=l 

If  there  are  exactly  enough  units  in  supply  to  meet  the  total  demand, 

m  n 

i=l  3=1 

then  every  solution  will  satisfy  the  inequalities  (C.2)  and  (C.3).  So,  to  solve  the 
problem  an  additional  demand  node,  n  +  1,  may  be  added  with  a  demand  equal  to 
the  excess  supply, 

m  n 

bn+l  ~  ^  ^  ^ 

i=\  3=1 

that  costs  nothing  to  deliver  i.e.  =  0. 

The  solution  approach  is  to  transform  the  problem  into  finding  the  maximum  flow 
through  a  network. 

C.2  Maximum  Flow  Through  a  Network 

Given  a  network,  N  =  {V,E),  there  is  a  source  node,  s,  and  a  sink  node,  t.  Each 
edge  has  a  capacity,  kij.  The  flow  through  any  edge,  Xij,  must  satisfy  0  <  Xij  <  kij. 
Furthermore,  all  the  nodes  other  than  s  and  t  must  obey  a  conservation  of  flow  given 
by 

^i,3  ~  X/  ~  b) 

i  I 

where  i  are  the  edges  flowing  into  node  j  and  I  are  the  edges  flowing  out  of  node  j. 
If  we  define  f{t)  as  the  flow  into  the  sink  node  t,  then  we  want  to  find  the  set  of  Xij 
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such  that  f{t)  is  a  maximum. 

To  find  the  maximum  fiow  in  a  network,  we  look  for  any  path  from  s  to  f  that 
can  augment  the  flow.  As  we  search  possible  edges  for  inclusion  into  an  augmenting 
path,  if  the  {i,  j?')th  edge  is  in  the  direction  of  the  flow,  then  Xij  must  be  less  than  kij 
for  this  edge  to  be  included;  if  the  edge  is  in  the  direction  opposing  the  flow  then  xij 
must  be  nonzero  to  be  included. 

To  find  a  path  from  s  to  t  that  augments  the  flow,  a  systematic  process  of  labeling 
and  scanning  the  nodes,  starting  at  s  and  ending  at  t,  is  performed.  At  any  give  time 
in  the  search,  a  node  can  be  in  one  of  three  states:  1.)  unlabeled  and  unscanned;  2.) 
labeled  and  unscanned;  or  3.)  labeled  and  scanned. 

Each  edge  has  a  capacity,  kij,  and  a  flow  Xij.  The  capacities  are  fixed  and  the  flows 
are  varied  to  achieve  the  maximum  flow.  All  flows  are  initially  set  to  zero  {xij  =  0). 
Thus,  in  the  beginning,  any  path  with  positive  flow  from  s  to  t  will  be  an  augmenting 
path.  The  labeling  of  a  node,  j,  assigns  two  numbers  to  it:  1.)  a  source  node,  i,  (not 
necessarily  s  and  2.)  the  maximum  augmentation  flow  possible,  fj,  from  that  source 
node. 

All  nodes  start  as  unlabeled  and  unscanned.  We  start  by  labeling  s  with  (— ,  oo). 
There  is  no  source  node  feeding  s  since  it  is  the  first  node  in  the  path  and  it  has 
an  infinite  capacity  to  augment  the  flow.  Node  s  is  now  considered  labeled  and 
unscanned. 

Next,  for  any  node  j  that  is  labeled  with  (i,  fj{or  -  fj))  and  unscanned,  the  node 
is  scanned  by  labeling  all  of  its  adjacent  nodes,  I,  that  are  currently  unlabeled.  If 
an  edge  is  directed  from  node  j  to  node  I  and  the  edge  has  excess  capacity  (ie. 
^i,j  <  h,j,)  then  we  label  node  I  with  {j,  fi)  where  f  is  the  maximum  amount  of  the 
flow  that  node  j  can  augment  to  node  I  through  edge  (j,  1).  The  amount  that  the 
flow  can  be  augmented  is  the  minimum  between  fj  (the  amount  at  node  j  available 
to  augment)  and  kij  —  Xij  (the  excess  capacity  available  to  augment). 

If  the  edge  is  directed  from  node  I  to  node  j  and  the  edge  has  nonzero  flow  {i.e. 
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Xij  >  0)  then  we  label  node  I  with  (j,  -fi)  where  fi  is  the  maximum  amount  of  the 
flow  that  node  j  can  augment  to  node  I  through  edge  {j,  1)  which  is  the  minimum 
between  fj  (the  amount  at  node  j  it  can  take  from  1)  and  Xij  (the  most  the  flow  can 
be  reduced  through  the  edge). 

The  labeling  and  scanning  continues  until  the  t  node  is  labeled  with  {I,  ft)  or  no 
path  can  be  found  to  t  to  augment  the  flow.  If  no  path  is  found,  the  flow  is  maximized 
and  the  problem  is  solved.  If  an  augmenting  path  is  found,  then  the  edge  flows  {xij) 
through  the  path  are  reassigned  to  increase  the  flow  from  s  to  thy  ft-  Next  all  the 
labels  are  removed  and  the  process  starts  over. 

There  is  a  theorem  in  Graph  Theory  called  Ford- Fulkerson  that  states  when  the 
labeling  algorithm  terminates  the  flow  is  optimal.  Furthermore,  if  the  algorithm 
always  uses  the  shortest  augmenting  path  {i.e.  fewest  edges)  the  time  to  solve  the 
problem  is  on  the  order  of  nrrF  where  n  is  the  number  of  nodes  and  m  is  the  number 
of  edges. 


C.3  Transportation  and  the  Maximum  Flow  Prob¬ 
lem 

The  transportation  problem  can  be  transformed  into  a  maximum  flow  problem  so 
that  the  algorithm  above  can  be  used. 

First,  we  insert  a  source  node  s  to  the  left  of  the  m  supply  nodes  Sj,  i  — 
1, 2, . . . ,  m.  A  directed  edge  with  capacity  at  is  drawn  from  s  to  Sj.  Next,  we  insert  a 
sink  node  t  to  the  right  of  the  n  demand  nodes  tj,  j  =  1,2, . . .  ,n.  A  directed  edge 
with  capacity  bi  is  drawn  from  L  to  t.  Figure  C-2  is  a  representation  of  the  network 
used  to  solve  the  transportation  problem.  Unlike  the  original  graph  not  all  edges 
between  the  source  and  demand  nodes  are  allowed. 

By  duality  in  Graph  Theory,  we  may  assign  dual  variables  Ui  and  Vj  to  the  con¬ 
straints  of  the  transportation  problem,  and  where  the  old  problem  minimizes  a  quan- 
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Figure  C-2:  Transportation  as  a  maximum  flow  network. 


tity,  the  new  problem  maximizes  a  new  quantity  based  on  the  old  constraints.  In  this 
case,  the  quantity 

m  n 

-  ^  Uitti  +  ^  Ujbj 
i=i  j=i 

is  maximized  within  the  constraints 


Ui  “f“  Vj  ^ 


Now  the  edges  from  Sj  to  tj  are  allowed  if  the  dual  variables  satisfy  —  Wj  +  Vj  =  Cij. 
Each  of  these  edges  has  infinite  capacity.  Given  a  set  of  Ui  and  Vj,  a  network  is 
constructed  and  the  maximum  flow  is  found.  However,  the  original  transportation 
problem  demand  constraints  must  be  checked  to  see  if  they  are  satisfied  i.e. 
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The  original  supply  constraints  are  inherently  satisfied  by  setting  the  capacity  between 
the  source  node  s  and  the  supply  node  Si  to  Oj. 

If  the  demand  constraint  is  met  the  problem  is  done;  otherwise,  a  new  set  of  Ui 
and  Vj  is  found  that  does  not  change  the  current  edge  flows.  In  other  words  new  edges 
must  be  established.  To  create  new  edges  we  define  a  flow  increment,  5, 

6  —  min[cj j  +  Ui  —  Vj  :  i  E  I,  jE  J] 

where  I  is  the  set  of  labeled  supply  nodes  and  J  is  the  set  of  unlabeled  demand  nodes. 
In  this  step  we’re  looking  for  the  lowest  cost  path  that  is  not  currently  being  used. 
Then  the  new  dual  variables  are  Ui  =  ui  if  node  Sj  is  labeled  and  ttj  =  Uj  +  6  if  Si  is 
unlabeled,  and  Vj  =  Vj  in  node  tj  is  labeled  and  Vj  =  Vj  +  S  if  tj  is  unlabeled. 

Then  the  maximum  flow  algorithm  is  run  and  the  demand  constraints  are  exam¬ 
ined  again.  The  process  repeats  until  the  demand  constraints  are  met. 


C.4  Application  to  Branch  Cut  Connections 

The  connection  of  positive  and  negative  residues  to  form  branch  cuts  is  a  very  special 
case  of  the  transportation  problem.  We  can  consider  the  positive  residues  as  the 
supply  nodes,  i,  and  the  negative  residues  as  the  demand  nodes,  j.  Since  there  are 
an  equal  number  of  positive  and  negative  residues,  m  =  n.  Furthermore,  since  a 
single  positive  residue  can  connect  to  exactly  one  negative  residue,  we  have  Oj  =  1 
and  bj  =  1.  Finally,  the  cost,  Cij,  of  transporting  one  unit  from  supply  node  i  to 
demand  node  j,  is  the  distance  (by  pixel  location)  between  the  corresponding  residue 
pair  defined  as  lij. 

Solving  the  transportation  problem  with 

di  —  1,  i  1,2,...,  777-, 
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is  equivalent  to  finding  the  set  of  branch  cut  connections  that  minimizes  the  sum  of 
all  branch  cut  lengths. 
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